module Statistics.MiniNest(
NestedSamplingResult(nsLogZ, nsLogZdelta, nsInfoNats, nsSamples),
SamplingObject(setLogWt, getLogWt, getLogL),
nestedSampling) where
import Control.Monad (forM)
import Data.IORef
import Data.List (sort)
import System.Random (randomRIO)
import Text.Printf
plus :: Double -> Double -> Double
plus x y
| x > y = x + log (1 + exp (yx))
| otherwise = y + log (1 + exp (xy))
data NestedSamplingResult a = NestedSamplingResult {
nsLogZ :: Double,
nsLogZdelta :: Double,
nsInfoNats :: Double,
nsSamples :: [a] }
instance Show (NestedSamplingResult a) where
show result =
(printf "logZ: %.2f +- %.2f\n" (nsLogZ result) (nsLogZdelta result) ++
printf "information: %.2f nats\n" (nsInfoNats result) ++
printf "%i samples\n" (length $ nsSamples result))
class SamplingObject a where
setLogWt :: a -> Double -> a
getLogWt :: a -> Double
getLogL :: a -> Double
nestedSampling :: (Ord a, SamplingObject a) => [a] -> (a -> Double -> IO a) -> Int -> IO (NestedSamplingResult a)
nestedSampling priorSamples explore iterations = do
let n = length priorSamples
objsRef <- newIORef priorSamples
samplesRef <- newIORef []
hRef <- newIORef 0
logZRef <- newIORef (10**37)
logWidthRef <- newIORef $ getLogWidth n
forM [1..iterations] (\nest -> do
objs <- readIORef objsRef
let worst = head $ sort objs
logwidth <- readIORef logWidthRef
let worst' = setLogWt worst (logwidth + (getLogL worst))
logZ <- readIORef logZRef
h <- readIORef hRef
let logZnew = plus logZ (getLogWt worst')
writeIORef hRef $ (exp $ getLogWt worst' logZnew) * (getLogL worst')
+ (exp $ logZ logZnew) * (h + logZ) logZnew
writeIORef logZRef logZnew
oldSamples <- readIORef samplesRef
writeIORef samplesRef (worst' : oldSamples)
let objs' = drop 1 $ sort objs
writeIORef objsRef objs'
objToCopy <- choice objs'
let logLstar = getLogL worst'
mutatedCopy <- explore objToCopy logLstar
writeIORef objsRef (mutatedCopy : objs')
writeIORef logWidthRef (logwidth 1.0 / fromIntegral n))
logZ <- readIORef logZRef
h <- readIORef hRef
samples <- readIORef samplesRef
return $ NestedSamplingResult {
nsLogZ=logZ,
nsLogZdelta=sqrt (h / fromIntegral n),
nsInfoNats=h,
nsSamples=samples
}
choice :: [a] -> IO a
choice [] = error "No items specified for choice."
choice [x] = return x
choice xs = do
let n = length xs
k <- randomRIO (0, n1)
return $ xs !! k
floatRatio :: Int -> Int -> Float
floatRatio n1 n2 = fromIntegral n1 / fromIntegral n2
getLogWidth :: Int -> Double
getLogWidth n = log $ 1.0 exp(1.0 / fromIntegral n)