Metropolis.hs (1238B)
1 2 module Deanie.Inference.Metropolis ( 3 metropolis 4 ) where 5 6 import Deanie.Language 7 import qualified Control.Foldl as L 8 import Control.Monad.Loops 9 10 data MHP a = MHP { 11 n :: {-# UNPACK #-} !Int 12 , current :: !a 13 , ccost :: {-# UNPACK #-} !Double 14 } 15 16 metropolis 17 :: Foldable f 18 => Int -> f a -> Program b -> (b -> a -> Double) 19 -> Program [b] 20 metropolis epochs obs prior model = do 21 current <- prior 22 unfoldrM mh MHP { n = epochs, ccost = cost current, .. } 23 where 24 cost param = L.fold (L.premap (model param) L.sum) obs 25 26 mh MHP {..} 27 | n <= 0 = return Nothing 28 | otherwise = do 29 proposal <- prior 30 let pcost = cost proposal 31 prob = moveProbability ccost pcost 32 accept <- bernoulli prob 33 let (nepochs, nlocation, ncost) = 34 if accept 35 then (pred n, proposal, pcost) 36 else (pred n, current, ccost) 37 return (Just (nlocation, MHP nepochs nlocation ncost)) 38 39 moveProbability :: Double -> Double -> Double 40 moveProbability current proposal = 41 whenNaN 0 (exp (min 0 (proposal - current))) 42 where 43 whenNaN val x 44 | isNaN x = val 45 | otherwise = x 46 {-# INLINE moveProbability #-} 47