commit 920c4203bb1742b894830e07003702cc00119f8c
parent 2c72f387f130d4db53d8316f5bfea52aeb3d97af
Author: Jared Tobin <jared@jtobin.ca>
Date: Mon, 26 Nov 2012 14:26:57 +1300
Slight API change. Update examples accordingly.
Diffstat:
7 files changed, 71 insertions(+), 115 deletions(-)
diff --git a/Examples/code/BNN_Flat.hs b/Examples/code/BNN_Flat.hs
@@ -23,19 +23,15 @@ main = do
inits <- readInits (args !! 1)
- let nepochs = read (head args) :: Int
- params = Options target (V.length inits) 30
- config = MarkovChain inits 0
-
- g <- create
- results <- runChain params nepochs 0 1 config g
-
- hPutStrLn stderr $
- let nAcc = accepts results
- total = nepochs * V.length inits * length (V.head inits)
- in show nAcc ++ " / " ++ show total ++ " (" ++
- show ((fromIntegral nAcc / fromIntegral total) :: Float) ++
- ") proposals accepted"
+ let nepochs = read (head args) :: Int
+ opts = Options { _size = V.length inits
+ , _nEpochs = nepochs
+ , _burnIn = 0
+ , _thinEvery = 1
+ , _csize = 30 }
+ initState = MarkovChain inits 0
+ g <- create
+ void $ runChain target opts initState g
diff --git a/Examples/code/Booth_Flat.hs b/Examples/code/Booth_Flat.hs
@@ -1,41 +0,0 @@
-import System.IO
-import System.Exit
-import System.Environment
-import System.Random.MWC
-import Control.Monad
-import Numeric.MCMC.Flat
-import qualified Data.Vector as V
-
-target :: [Double] -> Double
-target [x0, x1] = (-1)*((x0 + 2*x1 -7)^2 + (2*x0 + x1 - 5)^2)
-{-# INLINE target #-}
-
-main = do
- args <- getArgs
- when (args == []) $ do
- putStrLn "(flat-mcmc) Booth density "
- putStrLn "Usage: ./Booth_Flat <numSteps> <inits> "
- putStrLn " "
- putStrLn "numSteps : Number of Markov chain iterations to run."
- putStrLn "inits : Filepath containing points at which to "
- putStrLn " initialize the ensemble. "
- exitSuccess
-
- inits <- readInits (args !! 1)
-
- let nepochs = read (head args) :: Int
- params = Options target (V.length inits) 30
- config = MarkovChain inits 0
-
- g <- create
- results <- runChain params nepochs 0 1 config g
-
- hPutStrLn stderr $
- let nAcc = accepts results
- total = nepochs * V.length inits * length (V.head inits)
- in show nAcc ++ " / " ++ show total ++ " (" ++
- show ((fromIntegral nAcc / fromIntegral total) :: Float) ++
- ") proposals accepted"
-
-
-
diff --git a/Examples/code/Himmelblau_Flat.hs b/Examples/code/Himmelblau_Flat.hs
@@ -26,18 +26,15 @@ main = do
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 burnIn 1 config g
-
- hPutStrLn stderr $
- let nAcc = accepts results
- total = nepochs * V.length inits * length (V.head inits)
- in show nAcc ++ " / " ++ show total ++ " (" ++
- show ((fromIntegral nAcc / fromIntegral total) :: Float) ++
- ") proposals accepted"
+
+ opts = Options { _size = V.length inits
+ , _nEpochs = nepochs
+ , _burnIn = burnIn
+ , _thinEvery = 1
+ , _csize = 30 }
+ initState = MarkovChain inits 0
+ g <- create
+ void $ runChain target opts initState g
diff --git a/Examples/code/Rosenbrock_Flat.hs b/Examples/code/Rosenbrock_Flat.hs
@@ -24,16 +24,14 @@ main = do
inits <- readInits (args !! 1)
let nepochs = read (head args) :: Int
- params = Options target (V.length inits) 20
- config = MarkovChain inits 0
+ opts = Options { _size = V.length inits
+ , _nEpochs = nepochs
+ , _burnIn = 0
+ , _thinEvery = 1
+ , _csize = 30 }
- g <- create
- results <- runChain params nepochs 0 1 config g
+ initState = MarkovChain inits 0
- hPutStrLn stderr $
- let nAcc = accepts results
- total = nepochs * V.length inits * length (V.head inits)
- in show nAcc ++ " / " ++ show total ++ " (" ++
- show ((fromIntegral nAcc / fromIntegral total) :: Float) ++
- ") proposals accepted"
+ g <- create
+ void $ runChain target opts initState g
diff --git a/Examples/code/SPDE_Flat.hs b/Examples/code/SPDE_Flat.hs
@@ -34,16 +34,15 @@ main = do
thinEvery = read (args !! 2) :: Int
burnIn = read (args !! 3) :: Int
gran = read (args !! 4) :: Int
- params = Options target (V.length inits) gran
- config = MarkovChain inits 0
+
+ opts = Options { _size = V.length inits
+ , _nEpochs = nepochs
+ , _burnIn = burnIn
+ , _thinEvery = thinEvery
+ , _csize = gran }
+
+ initState = MarkovChain inits 0
g <- create
- results <- runChain params nepochs burnIn thinEvery config g
-
- hPutStrLn stderr $
- let nAcc = accepts results
- total = nepochs * V.length inits * length (V.head inits)
- in show nAcc ++ " / " ++ show total ++ " (" ++
- show ((fromIntegral nAcc / fromIntegral total) :: Float) ++
- ") proposals accepted"
+ void $ runChain target opts initState g
diff --git a/Numeric/MCMC/Flat.hs b/Numeric/MCMC/Flat.hs
@@ -40,9 +40,11 @@ instance Show MarkovChain where
-- | Options for the chain. The target (expected to be a log density), as
-- well as the size of the ensemble. The size should be an even number. Also
-- holds the specified parallel granularity as 'csize'.
-data Options = Options { _target :: [Double] -> Double
- , _size :: {-# UNPACK #-} !Int
- , _csize :: {-# UNPACK #-} !Int }
+data Options = Options { _size :: {-# UNPACK #-} !Int
+ , _nEpochs :: {-# UNPACK #-} !Int
+ , _burnIn :: {-# UNPACK #-} !Int
+ , _thinEvery :: {-# UNPACK #-} !Int
+ , _csize :: {-# UNPACK #-} !Int }
-- | An ensemble of particles.
type Ensemble = V.Vector [Double]
@@ -76,13 +78,14 @@ metropolisResult w0 w1 z zc target =
-- perturbing them with affine transformations based on particles
-- in a complementary ensemble, in parallel.
executeMoves :: (Functor m, PrimMonad m)
- => Ensemble -- Target sub-ensemble
+ => ([Double] -> Double) -- Target to sample
+ -> Ensemble -- Target sub-ensemble
-> Ensemble -- Complementary sub-ensemble
-> Int -- Size of the sub-ensembles
-> Gen (PrimState m) -- MWC PRNG
-> ViewsOptions m (Ensemble, Int) -- Updated ensemble and # of accepts
-executeMoves e0 e1 n g = do
- Options t _ csize <- ask
+executeMoves t e0 e1 n g = do
+ Options _ _ _ _ csize <- ask
zs <- replicateM n (lift $ symmetricVariate g)
zcs <- replicateM n (lift $ uniformR (0 :: Double, 1 :: Double) g)
@@ -103,18 +106,19 @@ executeMoves e0 e1 n g = do
-- perturbing each element and accepting/rejecting the perturbation in
-- parallel.
metropolisStep :: (Functor m, PrimMonad m)
- => MarkovChain -- State of the Markov chain
+ => ([Double] -> Double) -- Target to sample
+ -> MarkovChain -- State of the Markov chain
-> Gen (PrimState m) -- MWC PRNG
-> ViewsOptions m MarkovChain -- Updated sub-ensemble
-metropolisStep state g = do
- Options _ n _ <- ask
+metropolisStep t state g = do
+ Options n _ _ _ _ <- ask
let n0 = truncate (fromIntegral n / (2 :: Double)) :: Int
(e, nacc) = (ensemble &&& accepts) state
(e0, e1) = (V.slice (0 :: Int) n0 &&& V.slice n0 n0) e
-- Update each sub-ensemble
- result0 <- executeMoves e0 e1 n0 g
- result1 <- executeMoves e1 (fst result0) n0 g
+ result0 <- executeMoves t e0 e1 n0 g
+ result1 <- executeMoves t e1 (fst result0) n0 g
return $!
MarkovChain (V.concat $ map fst [result0, result1])
@@ -122,33 +126,36 @@ metropolisStep state g = do
{-# INLINE metropolisStep #-}
-- | 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 burnIn thinEvery initConfig g
+runChain :: ([Double] -> Double) -- ^ Target to sample
+ -> Options -- ^ Options of the Markov chain
+ -> MarkovChain -- ^ Initial state of the Markov chain
+ -> Gen RealWorld -- ^ MWC PRNG
+ -> IO MarkovChain -- ^ End state of the Markov chain, wrapped in IO
+runChain target opts initState g
| l == 0
= error "runChain: ensemble must contain at least one particle"
- | l < (length . V.head) (ensemble initConfig)
+ | l < (length . V.head) (ensemble initState)
= do hPutStrLn stderr $ "runChain: ensemble should be twice as large as "
++ "the target's dimension. Continuing anyway."
- go opts nepochs thinEvery initConfig g
+ go opts nepochs thinEvery initState g
| burnIn < 0 || thinEvery < 0 = error "runChain: nonsensical burn-in or thinning input."
- | otherwise = go opts nepochs thinEvery initConfig g
+ | otherwise = go opts nepochs thinEvery initState g
where
- l = V.length (ensemble initConfig)
- go o n t !c g0 | n == 0 = return c
+ Options l nepochs burnIn thinEvery _ = opts
+ go o n t !c g0 | n == 0 = hPutStrLn stderr
+ (let nAcc = accepts c
+ total = nepochs * l * length (V.head $ ensemble c)
+ in show nAcc ++ " / " ++ show total ++ " (" ++
+ show ((fromIntegral nAcc / fromIntegral total) :: Float) ++
+ ") proposals accepted") >> return c
| n > (nepochs - burnIn) = do
- r <- runReaderT (metropolisStep c g0) o
+ r <- runReaderT (metropolisStep target c g0) o
go o (n - 1) t r g0
| n `rem` t /= 0 = do
- r <- runReaderT (metropolisStep c g0) o
+ r <- runReaderT (metropolisStep target c g0) o
go o (n - 1) t r g0
| otherwise = do
- r <- runReaderT (metropolisStep c g0) o
+ r <- runReaderT (metropolisStep target c g0) o
print r
go o (n - 1) t r g0
{-# INLINE runChain #-}
diff --git a/flat-mcmc.cabal b/flat-mcmc.cabal
@@ -2,7 +2,7 @@
-- documentation, see http://haskell.org/cabal/users-guide/
name: flat-mcmc
-version: 0.1.0.0
+version: 0.2.0.0
synopsis: Painless general-purpose sampling.
-- description:
homepage: http://jtobin.github.com/flat-mcmc
@@ -26,6 +26,6 @@ Source-repository head
library
exposed-modules: Numeric.MCMC.Flat
-- other-modules:
- build-depends: base ==4.*, mtl ==2.1.*, primitive ==0.4.*, mwc-random ==0.12.*, vector ==0.9.*, monad-par ==0.3.*, monad-par-extras ==0.3.*
+ build-depends: base >= 4.3, mtl >= 2.1, primitive >= 0.4, mwc-random >= 0.12, vector >= 0.9, monad-par >= 0.3, monad-par-extras >= 0.3
ghc-options: -Wall
