multiple chains in parallel

[aezarebski]

Running multiple chains in parallel is a good way to reduce Monte Carlo error in estimates and provides simple diagnostics for convergence and mixing. However it is non-trivial to run multiple chains in parallel. A function multipleChains which took care of the details of running multiple chains in parallel would be nice. Perhaps pipes-concurrency would be useful in implementing this.

[aezarebski]

Another approach to this is to use the functionality from monad-par. By generating the samples in the ST monad we can avoid the IO monad which seems to simplify things. Here is an example that should run two chains in parallel.

module Main where

import Numeric.MCMC.Metropolis
import qualified System.Random.MWC as MWC
import Control.Monad.ST
import Control.Monad.Par (runPar,Par(..))
import Control.Monad.Par.Combinator (parMap)

rosenbrock :: [Double] -> Double
rosenbrock [x0, x1] = negate (5  *(x1 - x0 ** 2.0) ** 2.0 + 0.05 * (1 - x0) ** 2.0)


singleChain :: [Double] -> [Chain [Double] b0]
singleChain ix = let f = MWC.asGenST $ \gen -> chain 10 1 ix rosenbrock gen
         in runST $ do genny <- MWC.create
                       f genny

parallelExpression :: Par [[[Double]]]
parallelExpression = parMap (drop 5 . map chainPosition . singleChain) [[0,0],[1,1]]

main :: IO ()
main = do
  putStrLn "Hello"
  print $ runPar parallelExpression
  putStrLn "Goodbye"

this does require a couple of packages to be added: monad-par, monad-par-extras and mwc-random. A problem with this approach is that each chain uses the same seed.