flat-mcmc

commit 2c72f387f130d4db53d8316f5bfea52aeb3d97af
parent 3497d940f55edd5e6e745dd96559b0b440b9a82e
Author: Jared Tobin <jared@jtobin.ca>
Date:   Mon, 12 Nov 2012 10:33:30 +1300

Add burn-in option to runChain.

Diffstat:
M
Examples/code/BNN_Flat.hs
 | 
2
+-
M
Examples/code/Booth_Flat.hs
 | 
2
+-
M
Examples/code/Himmelblau_Flat.hs
 | 
8
+++++---
M
Examples/code/Rosenbrock_Flat.hs
 | 
2
+-
M
Examples/code/SPDE_Flat.hs
 | 
10
+++++-----
M
Numeric/MCMC/Flat.hs
 | 
7
++++++-

6 files changed, 19 insertions(+), 12 deletions(-)
diff --git a/Examples/code/BNN_Flat.hs b/Examples/code/BNN_Flat.hs
@@ -28,7 +28,7 @@ main = do
         config  = MarkovChain inits 0
 
     g       <- create
-    results <- runChain params nepochs 1 config g
+    results <- runChain params nepochs 0 1 config g
 
     hPutStrLn stderr $ 
         let nAcc  = accepts results
diff --git a/Examples/code/Booth_Flat.hs b/Examples/code/Booth_Flat.hs
@@ -28,7 +28,7 @@ main = do
         config  = MarkovChain inits 0
 
     g       <- create
-    results <- runChain params nepochs 1 config g
+    results <- runChain params nepochs 0 1 config g
 
     hPutStrLn stderr $ 
         let nAcc  = accepts results
diff --git a/Examples/code/Himmelblau_Flat.hs b/Examples/code/Himmelblau_Flat.hs
@@ -14,21 +14,23 @@ main = do
     args  <- getArgs 
     when (args == []) $ do
         putStrLn  "(flat-mcmc) Himmelblau density                              "
-        putStrLn  "Usage: ./Himmelblau_Flat <numSteps> <inits>                 " 
+        putStrLn  "Usage: ./Himmelblau_Flat <numSteps> <burnIn> <inits>        " 
         putStrLn  "                                                            "
         putStrLn  "numSteps         : Number of Markov chain iterations to run."
+        putStrLn  "burnIn           : Number of burn-in steps to perform.      "
         putStrLn  "inits            : Filepath containing points at which to   "
         putStrLn  "                   initialize the ensemble.                 "
         exitSuccess
 
-    inits <- readInits (args !! 1)
+    inits <- readInits (args !! 2)
 
     let nepochs = read (head args) :: Int
+        burnIn  = read (args !! 1) :: Int
         params  = Options     target (V.length inits) 30
         config  = MarkovChain inits 0
 
     g       <- create
-    results <- runChain params nepochs 1 config g
+    results <- runChain params nepochs burnIn 1 config g
 
     hPutStrLn stderr $ 
         let nAcc  = accepts results
diff --git a/Examples/code/Rosenbrock_Flat.hs b/Examples/code/Rosenbrock_Flat.hs
@@ -28,7 +28,7 @@ main = do
         config = MarkovChain inits 0
 
     g       <- create
-    results <- runChain params nepochs 1 config g
+    results <- runChain params nepochs 0 1 config g
 
     hPutStrLn stderr $ 
         let nAcc  = accepts results
diff --git a/Examples/code/SPDE_Flat.hs b/Examples/code/SPDE_Flat.hs
@@ -1,5 +1,3 @@
-{-# LANGUAGE BangPatterns #-}
-
 import System.IO
 import System.Exit
 import System.Environment
@@ -20,12 +18,13 @@ main = do
     args  <- getArgs 
     when (args == []) $ do
         putStrLn  "(flat-mcmc) Stochastic partial differential equation        "
-        putStrLn  "Usage: ./SPDE_Flat <numSteps> <inits> <thinEvery> <granularity>         " 
+        putStrLn  "Usage: ./SPDE_Flat <numSteps> <inits> <thinEvery> <burnIn> <granularity>         " 
         putStrLn  "                                                            "
         putStrLn  "numSteps         : Number of Markov chain iterations to run."
         putStrLn  "inits            : Filepath containing points at which to   "
         putStrLn  "                   initialize the ensemble.                 "
         putStrLn  "thinEvery        : Print every n^th iteration.              "
+        putStrLn  "burnIn           : Number of burn-in steps.                 "
         putStrLn  "granularity      : Parallel granularity (smaller is finer). "
         exitSuccess
 
@@ -33,12 +32,13 @@ main = do
 
     let nepochs   = read (head args) :: Int
         thinEvery = read (args !! 2) :: Int
-        gran      = read (args !! 3) :: Int
+        burnIn    = read (args !! 3) :: Int
+        gran      = read (args !! 4) :: Int
         params    = Options target (V.length inits) gran
         config    = MarkovChain inits 0
 
     g       <- create
-    results <- runChain params nepochs thinEvery config g
+    results <- runChain params nepochs burnIn thinEvery config g
 
     hPutStrLn stderr $ 
         let nAcc  = accepts results
diff --git a/Numeric/MCMC/Flat.hs b/Numeric/MCMC/Flat.hs
@@ -124,21 +124,26 @@ metropolisStep state g = do
 -- | Diffuse through states.
 runChain :: Options         -- Options of the Markov chain
          -> Int             -- Number of epochs to iterate the chain
+         -> Int             -- Burn-in period
          -> 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 opts nepochs thinEvery initConfig g 
+runChain opts nepochs burnIn thinEvery initConfig g 
     | l == 0 
         = error "runChain: ensemble must contain at least one particle"
     | l < (length . V.head) (ensemble initConfig)
         = do hPutStrLn stderr $ "runChain: ensemble should be twice as large as "
                              ++ "the target's dimension.  Continuing anyway."
              go opts nepochs thinEvery initConfig g
+    | burnIn < 0 || thinEvery < 0 = error "runChain: nonsensical burn-in or thinning input."
     | otherwise = go opts nepochs thinEvery initConfig g
   where 
     l = V.length (ensemble initConfig)
     go o n t !c g0 | n == 0 = return c
+                   | n > (nepochs - burnIn) = do
+                       r <- runReaderT (metropolisStep c g0) o
+                       go o (n - 1) t r g0 
                    | n `rem` t /= 0 = do
                        r <- runReaderT (metropolisStep c g0) o
                        go o (n - 1) t r g0