diff --git a/NestedSampling.cabal b/NestedSampling.cabal --- a/NestedSampling.cabal +++ b/NestedSampling.cabal @@ -1,5 +1,5 @@ Name: NestedSampling -Version: 0.1.2 +Version: 0.1.3 Synopsis: A port of John Skilling's nested sampling C code to Haskell. Description: Nested Sampling is a numerical algorithm for approximate Bayesian @@ -29,15 +29,14 @@ Copyright: (C) Sivia, Skilling 2006, Trotts 2011 Category: Statistics Build-type: Simple -Extra-source-files: lighthouse.hs README -Cabal-version: >=1.8 +Extra-source-files: examples/lighthouse.hs README +Cabal-version: >=1.6 +source-repository head + type: git + location: git://github.com/ijt/haskell_nested_sampling.git Library Exposed-modules: Statistics.MiniNest - Build-depends: base >= 4 && < 5, random + Build-depends: base >= 4 && < 5, random, vector hs-source-dirs: lib - -Executable lighthouse - Main-Is: lighthouse.hs - Build-depends: base >= 4 && < 5, random, vector, NestedSampling diff --git a/examples/lighthouse.hs b/examples/lighthouse.hs new file mode 100644 --- /dev/null +++ b/examples/lighthouse.hs @@ -0,0 +1,159 @@ +#!/usr/bin/env runhaskell + +-- lighthouse.hs "LIGHTHOUSE" NESTED SAMPLING APPLICATION +-- (GNU General Public License software, (C) Sivia and Skilling 2006) +-- u=0 u=1 +-- ------------------------------------- +-- y=2 |:::::::::::::::::::::::::::::::::::::| v=1 +-- |::::::::::::::::::::::LIGHT::::::::::| +-- north|::::::::::::::::::::::HOUSE::::::::::| +-- |:::::::::::::::::::::::::::::::::::::| +-- |:::::::::::::::::::::::::::::::::::::| +-- y=0 |:::::::::::::::::::::::::::::::::::::| v=0 +-- --*--------------*----*--------*-**--**--*-*-------------*-------- +-- x=-2 coastline -->east x=2 +-- Problem: +-- Lighthouse at (x,y) emitted n flashes observed at D[.] on coast. +-- Inputs: +-- Prior(u) is uniform (=1) over (0,1), mapped to x = 4*u - 2; and +-- Prior(v) is uniform (=1) over (0,1), mapped to y = 2*v; so that +-- Position is 2-dimensional -2 < x < 2, 0 < y < 2 with flat prior +-- Likelihood is L(x,y) = PRODUCT[k] (y/pi) / ((D[k] - x)^2 + y^2) +-- Outputs: +-- Evidence is Z = INTEGRAL L(x,y) Prior(x,y) dxdy +-- Posterior is P(x,y) = L(x,y) / Z estimating lighthouse position +-- Information is H = INTEGRAL P(x,y) log(P(x,y)/Prior(x,y)) dxdy + +import qualified Data.Vector.Unboxed as UV +import Control.Monad (mapM) +import Statistics.MiniNest +import System.Random (randomIO) +import Text.Printf + +data Lighthouse = Lighthouse { + lhU :: Double, + lhV :: Double, + lhX :: Double, + lhY :: Double, + lhLogL :: Double, + lhLogWt :: Double +} deriving (Eq, Show) + +instance Ord Lighthouse where + a <= b = lhLogL a <= lhLogL b + +instance SamplingObject Lighthouse where + setLogWt lh newLogWt = lh { lhLogWt = newLogWt } + getLogWt lh = lhLogWt lh + getLogL lh = lhLogL lh + +logLhoodOfData :: UV.Vector Double -> Double -> Double -> Double +logLhoodOfData observations x y = UV.sum $ UV.map term observations + where term dk = log (y / pi) - log ((dk - x)*(dk - x) + y*y) + +-- logLikelihood function +-- x: Easterly position +-- y: Northerly position +logLhood :: Double -> Double -> Double +logLhood x y = logLhoodOfData lhData x y + +lhData = UV.fromList [4.73, 0.45, -1.73, 1.09, 2.19, 0.12, 1.31, + 1.00, 1.32, 1.07, 0.86, -0.49, -2.59, 1.73, 2.11, + 1.61, 4.98, 1.71, 2.23,-57.20, 0.96, 1.25, -1.56, + 2.45, 1.19, 2.17,-10.66, 1.91, -4.16, 1.92, 0.10, 1.98, + -2.51, 5.55, -0.47, 1.91, 0.95, -0.78, -0.84, 1.72, + -0.01, 1.48, 2.70, 1.21, 4.41, -4.79, 1.33, 0.81, + 0.20, 1.58, 1.29, 16.19, 2.75, -2.38, -1.79, + 6.50,-18.53, 0.72, 0.94, 3.64, 1.94, -0.11, 1.57, 0.57] + +-- |Sample from U[0,1] +uniform :: IO Double +uniform = randomIO + +sampleFromPrior :: IO Lighthouse +sampleFromPrior = do + u <- uniform + v <- uniform + let x=4*u - 2 + y=2*v + return $ Lighthouse u v x y (logLhood x y) 0 + +-- |Evolve Lighthouse within likelihood constraint +-- obj: Lighthouse being evolved +-- logLstar: Likelihood constraint L > Lstar +explore :: Lighthouse -> Double -> IO Lighthouse +explore obj logLstar = + explore' step m accept reject (lhU obj) (lhV obj) (lhX obj) (lhY obj) + (lhLogL obj) + where step = 0.1 -- Initial guess suitable step-size in (0,1) + m = 20 -- MCMC counter (pre-judged # steps) + accept = 0 -- # MCMC acceptances + reject = 0 -- # MCMC rejections + explore' step m accept reject u v x y logL = do + -- Trial Lighthouse + unif1 <- uniform + unif2 <- uniform + let u' = wrapAround $ u + step * (2*unif1 - 1) -- |move| < step + v' = wrapAround $ v + step * (2*unif2 - 1) -- |move| < step + x' = 4*u' - 2 -- map to x + y' = 2*v' -- map to y + logL' = logLhood x' y' + + -- Accept if and only if within hard likelihood constraint + obj' <- + if logL' > logLstar + then return $ Lighthouse u' v' x' y' logL' (lhLogWt obj) + else return $ Lighthouse u v x y logL (lhLogWt obj) + (accept, reject) <- if logL' > logLstar + then return (accept + 1, reject) + else return (accept, reject + 1) + + -- Refine step-size to let acceptance ratio converge around 50% + step <- if accept > reject + then return $ step * exp(1.0 / accept) + else return step + step <- if accept < reject + then return $ step / exp(1.0 / reject) + else return step + if m == 0 + then return obj' + else explore' step (m-1) accept reject (lhU obj') (lhV obj') + (lhX obj') (lhY obj') (lhLogL obj') + +wrapAround :: Double -> Double +wrapAround x = x - (fromIntegral $ floor x) + +data Stats = Stats { meanX :: Double, + meanY :: Double, + stddevX :: Double, + stddevY :: Double } + +instance Show Stats where + show s = (printf "x = %.2f +- %.2f\n" (meanX s) (stddevX s) ++ + printf "y = %.2f +- %.2f\n" (meanY s) (stddevY s)) + +-- Posterior properties, here mean and stddev of x,y +-- Args: +-- samples: Objects defining posterior +-- logZ: Evidence (= total weight = SUM[Samples] Weight) +getStats :: [Lighthouse] -> Double -> Stats +getStats samples logZ = + Stats {meanX=x, + meanY=y, + stddevX=sqrt $ xx - x*x, + stddevY=sqrt $ yy - y*y } + where weightsSamples = [(exp (lhLogWt s - logZ), s) | s <- samples] + x = sum [w*(lhX s) | (w,s) <- weightsSamples] + y = sum [w*(lhY s) | (w,s) <- weightsSamples] + xx = sum [w*(lhX s)^2 | (w,s) <- weightsSamples] + yy = sum [w*(lhY s)^2 | (w,s) <- weightsSamples] + +main = do + let n = 100 -- # number of candidate lighthouses + let maxIterations = 1000 -- # iterates + priorSamples <- mapM (\_ -> sampleFromPrior) [1..n] + result <- nestedSampling priorSamples explore maxIterations + let stats = getStats (nsSamples result) (nsLogZ result) + print result + print stats + diff --git a/lighthouse.hs b/lighthouse.hs deleted file mode 100644 --- a/lighthouse.hs +++ /dev/null @@ -1,157 +0,0 @@ --- lighthouse.hs "LIGHTHOUSE" NESTED SAMPLING APPLICATION --- (GNU General Public License software, (C) Sivia and Skilling 2006) --- u=0 u=1 --- ------------------------------------- --- y=2 |:::::::::::::::::::::::::::::::::::::| v=1 --- |::::::::::::::::::::::LIGHT::::::::::| --- north|::::::::::::::::::::::HOUSE::::::::::| --- |:::::::::::::::::::::::::::::::::::::| --- |:::::::::::::::::::::::::::::::::::::| --- y=0 |:::::::::::::::::::::::::::::::::::::| v=0 --- --*--------------*----*--------*-**--**--*-*-------------*-------- --- x=-2 coastline -->east x=2 --- Problem: --- Lighthouse at (x,y) emitted n flashes observed at D[.] on coast. --- Inputs: --- Prior(u) is uniform (=1) over (0,1), mapped to x = 4*u - 2; and --- Prior(v) is uniform (=1) over (0,1), mapped to y = 2*v; so that --- Position is 2-dimensional -2 < x < 2, 0 < y < 2 with flat prior --- Likelihood is L(x,y) = PRODUCT[k] (y/pi) / ((D[k] - x)^2 + y^2) --- Outputs: --- Evidence is Z = INTEGRAL L(x,y) Prior(x,y) dxdy --- Posterior is P(x,y) = L(x,y) / Z estimating lighthouse position --- Information is H = INTEGRAL P(x,y) log(P(x,y)/Prior(x,y)) dxdy - -import qualified Data.Vector.Unboxed as UV -import Control.Monad (mapM) -import Statistics.MiniNest -import System.Random (randomIO) -import Text.Printf - -data Lighthouse = Lighthouse { - lhU :: Double, - lhV :: Double, - lhX :: Double, - lhY :: Double, - lhLogL :: Double, - lhLogWt :: Double -} deriving (Eq, Show) - -instance Ord Lighthouse where - a <= b = lhLogL a <= lhLogL b - -instance SamplingObject Lighthouse where - setLogWt lh newLogWt = lh { lhLogWt = newLogWt } - getLogWt lh = lhLogWt lh - getLogL lh = lhLogL lh - -logLhoodOfData :: UV.Vector Double -> Double -> Double -> Double -logLhoodOfData observations x y = UV.sum $ UV.map term observations - where term dk = log (y / pi) - log ((dk - x)*(dk - x) + y*y) - --- logLikelihood function --- x: Easterly position --- y: Northerly position -logLhood :: Double -> Double -> Double -logLhood x y = logLhoodOfData lhData x y - -lhData = UV.fromList [4.73, 0.45, -1.73, 1.09, 2.19, 0.12, 1.31, - 1.00, 1.32, 1.07, 0.86, -0.49, -2.59, 1.73, 2.11, - 1.61, 4.98, 1.71, 2.23,-57.20, 0.96, 1.25, -1.56, - 2.45, 1.19, 2.17,-10.66, 1.91, -4.16, 1.92, 0.10, 1.98, - -2.51, 5.55, -0.47, 1.91, 0.95, -0.78, -0.84, 1.72, - -0.01, 1.48, 2.70, 1.21, 4.41, -4.79, 1.33, 0.81, - 0.20, 1.58, 1.29, 16.19, 2.75, -2.38, -1.79, - 6.50,-18.53, 0.72, 0.94, 3.64, 1.94, -0.11, 1.57, 0.57] - --- |Sample from U[0,1] -uniform :: IO Double -uniform = randomIO - -sampleFromPrior :: IO Lighthouse -sampleFromPrior = do - u <- uniform - v <- uniform - let x=4*u - 2 - y=2*v - return $ Lighthouse u v x y (logLhood x y) 0 - --- |Evolve Lighthouse within likelihood constraint --- obj: Lighthouse being evolved --- logLstar: Likelihood constraint L > Lstar -explore :: Lighthouse -> Double -> IO Lighthouse -explore obj logLstar = - explore' step m accept reject (lhU obj) (lhV obj) (lhX obj) (lhY obj) - (lhLogL obj) - where step = 0.1 -- Initial guess suitable step-size in (0,1) - m = 20 -- MCMC counter (pre-judged # steps) - accept = 0 -- # MCMC acceptances - reject = 0 -- # MCMC rejections - explore' step m accept reject u v x y logL = do - -- Trial Lighthouse - unif1 <- uniform - unif2 <- uniform - let u' = wrapAround $ u + step * (2*unif1 - 1) -- |move| < step - v' = wrapAround $ v + step * (2*unif2 - 1) -- |move| < step - x' = 4*u' - 2 -- map to x - y' = 2*v' -- map to y - logL' = logLhood x' y' - - -- Accept if and only if within hard likelihood constraint - obj' <- - if logL' > logLstar - then return $ Lighthouse u' v' x' y' logL' (lhLogWt obj) - else return $ Lighthouse u v x y logL (lhLogWt obj) - (accept, reject) <- if logL' > logLstar - then return (accept + 1, reject) - else return (accept, reject + 1) - - -- Refine step-size to let acceptance ratio converge around 50% - step <- if accept > reject - then return $ step * exp(1.0 / accept) - else return step - step <- if accept < reject - then return $ step / exp(1.0 / reject) - else return step - if m == 0 - then return obj' - else explore' step (m-1) accept reject (lhU obj') (lhV obj') - (lhX obj') (lhY obj') (lhLogL obj') - -wrapAround :: Double -> Double -wrapAround x = x - (fromIntegral $ floor x) - -data Stats = Stats { meanX :: Double, - meanY :: Double, - stddevX :: Double, - stddevY :: Double } - -instance Show Stats where - show s = (printf "x = %.2f +- %.2f\n" (meanX s) (stddevX s) ++ - printf "y = %.2f +- %.2f\n" (meanY s) (stddevY s)) - --- Posterior properties, here mean and stddev of x,y --- Args: --- samples: Objects defining posterior --- logZ: Evidence (= total weight = SUM[Samples] Weight) -getStats :: [Lighthouse] -> Double -> Stats -getStats samples logZ = - Stats {meanX=x, - meanY=y, - stddevX=sqrt $ xx - x*x, - stddevY=sqrt $ yy - y*y } - where weightsSamples = [(exp (lhLogWt s - logZ), s) | s <- samples] - x = sum [w*(lhX s) | (w,s) <- weightsSamples] - y = sum [w*(lhY s) | (w,s) <- weightsSamples] - xx = sum [w*(lhX s)^2 | (w,s) <- weightsSamples] - yy = sum [w*(lhY s)^2 | (w,s) <- weightsSamples] - -main = do - let n = 100 -- # number of candidate lighthouses - let maxIterations = 1000 -- # iterates - priorSamples <- mapM (\_ -> sampleFromPrior) [1..n] - result <- nestedSampling priorSamples explore maxIterations - let stats = getStats (nsSamples result) (nsLogZ result) - print result - print stats -