NestedSampling 0.1.2 → 0.1.3
raw patch · 3 files changed
+166/−165 lines, 3 filesdep −NestedSamplingPVP ok
version bump matches the API change (PVP)
Dependencies removed: NestedSampling
API changes (from Hackage documentation)
Files
- NestedSampling.cabal +7/−8
- examples/lighthouse.hs +159/−0
- lighthouse.hs +0/−157
NestedSampling.cabal view
@@ -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
+ examples/lighthouse.hs view
@@ -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+
− lighthouse.hs
@@ -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-