commit 844e51dae3177d7b909eec7999f154322afe05e9
parent 22ea3f7d2baf4107774128ea3fad2fb0d0e173de
Author: Jared Tobin <jared@jtobin.ca>
Date: Mon, 5 Nov 2012 09:19:02 +1300
Add thinning parameter.
Diffstat:
2 files changed, 29 insertions(+), 21 deletions(-)
diff --git a/Examples/code/Rosenbrock_HMC.hs b/Examples/code/Rosenbrock_HMC.hs
@@ -15,26 +15,28 @@ gTarget = grad target
main = do
args <- getArgs
when (args == []) $ do
- putStrLn "(hasty-hamiltonian) Rosenbrock density "
- putStrLn "Usage: ./Rosenbrock_HMC <numSteps> <nDisc> <stepSize> <inits>"
- putStrLn " "
- putStrLn "numSteps : Number of Markov chain iterations to run. "
- putStrLn "nDisc : Number of discretizing steps to take. "
- putStrLn "stepSize : Perturbation scaling parameter. "
- putStrLn "inits : Filepath containing points at which to "
- putStrLn " initialize the chain. "
+ putStrLn "(hasty-hamiltonian) Rosenbrock density "
+ putStrLn "Usage: ./Rosenbrock_HMC <numSteps> <thinEvery> <nDisc> <stepSize> <inits>"
+ putStrLn " "
+ putStrLn "numSteps : Number of Markov chain iterations to run. "
+ putStrLn "thinEvery : Print every nth iteration. "
+ putStrLn "nDisc : Number of discretizing steps to take. "
+ putStrLn "stepSize : Perturbation scaling parameter. "
+ putStrLn "inits : Filepath containing points at which to "
+ putStrLn " initialize the chain. "
exitSuccess
inits <- fmap (map read . words) (readFile (args !! 3)) :: IO [Double]
- let nepochs = read (head args) :: Int
- nDisc = read (args !! 1) :: Double
- eps = read (args !! 2) :: Double
- params = Options target gTarget nDisc eps
- config = MarkovChain inits 0
+ let nepochs = read (head args) :: Int
+ thinEvery = read (args !! 1) :: Int
+ nDisc = read (args !! 2) :: Double
+ eps = read (args !! 3) :: Double
+ params = Options target gTarget nDisc eps
+ config = MarkovChain inits 0
g <- create
- results <- runChain params nepochs config g
+ results <- runChain params nepochs thinEvery config g
hPutStrLn stderr $
let nAcc = accepts results
diff --git a/Numeric/MCMC/Hamiltonian.hs b/Numeric/MCMC/Hamiltonian.hs
@@ -41,6 +41,7 @@ leapfrog t0 r0 = do
rt = zipWith (+) rm (map (* (0.5*e)) (gTarget tt))
in go tt rt (n - 1)
return $! go t0 r0 ndisc
+{-# INLINE leapfrog #-}
-- | Perform a Metropolis accept/reject step.
metropolisStep :: PrimMonad m
@@ -64,18 +65,23 @@ metropolisStep state g = do
- target t0 - 0.5*(r <.> r)
return $! MarkovChain (fst mc) (nacc + snd mc)
+{-# INLINE metropolisStep #-}
-- | Diffuse through states.
runChain :: Options -- Options of the Markov chain.
-> Int -- Number of epochs to iterate the chain.
+ -> Int -- Print every nth iteration.
-> MarkovChain -- Initial state of the Markov chain.
-> Gen RealWorld -- MWC PRNG
-> IO MarkovChain -- End state of the Markov chain, wrapped in IO.
-runChain params nepochs initConfig g
- | nepochs == 0 = return initConfig
- | otherwise = do
- result <- runReaderT (metropolisStep initConfig g) params
- print result
- runChain params (nepochs - 1) result g
+runChain = go
+ where go o n t !c g | n == 0 = return c
+ | n `rem` t /= 0 = do
+ r <- runReaderT (metropolisStep c g) o
+ go o (n - 1) t r g
+ | otherwise = do
+ r <- runReaderT (metropolisStep c g) o
+ print r
+ go o (n - 1) t r g
+{-# INLINE runChain #-}
-
