Copyright | (c) 2016 Jared Tobin |
---|---|
License | MIT |
Maintainer | Jared Tobin <jared@jtobin.ca> |
Stability | unstable |
Portability | ghc |
Safe Haskell | None |
Language | Haskell2010 |
This is the 'affine invariant ensemble' or AIEMCMC algorithm described in Goodman and Weare, 2010. It is a reasonably efficient and hassle-free sampler, requiring no mucking with tuning parameters or local proposal distributions.
The mcmc
function streams a trace to stdout to be processed elsewhere,
while the flat
transition can be used for more flexible purposes,
such as working with samples in memory.
See http://msp.org/camcos/2010/5-1/camcos-v5-n1-p04-p.pdf for the definitive reference of the algorithm.
- mcmc :: Int -> Ensemble -> (Particle -> Double) -> Gen RealWorld -> IO ()
- flat :: PrimMonad m => Transition m Chain
- type Particle = Vector Double
- type Ensemble = Vector Particle
- data Chain
- create :: PrimMonad m => m (Gen (PrimState m))
- createSystemRandom :: IO GenIO
- withSystemRandom :: PrimBase m => (Gen (PrimState m) -> m a) -> IO a
- asGenIO :: (GenIO -> IO a) -> GenIO -> IO a
Documentation
mcmc :: Int -> Ensemble -> (Particle -> Double) -> Gen RealWorld -> IO () Source
Trace n
iterations of a Markov chain and stream them to stdout.
Note that the Markov chain is defined over the space of ensembles, so you'll need to provide an ensemble of particles for the start location.
>>>
import Numeric.MCMC.Flat
>>>
import Data.Vector (Vector, toList, fromList)
>>>
:{
>>>
let rosenbrock xs = negate (5 *(x1 - x0 ^ 2) ^ 2 + 0.05 * (1 - x0) ^ 2)
where [x0, x1] = toList xs>>>
:}
>>>
:{
>>>
let ensemble = fromList [
>>>
fromList [negate 1.0, negate 1.0]
>>>
, fromList [negate 1.0, 1.0]
>>>
, fromList [1.0, negate 1.0]
>>>
, fromList [1.0, 1.0]
>>>
]
>>>
:}
>>>
withSystemRandom . asGenIO $ mcmc 2 ensemble rosenbrock
-1.0,-1.0 -1.0,1.0 1.0,-1.0 0.7049046915549257,0.7049046915549257 -0.843493377618159,-0.843493377618159 -1.1655594505975082,1.1655594505975082 0.5466534497342876,-0.9615123448709006 0.7049046915549257,0.7049046915549257
flat :: PrimMonad m => Transition m Chain Source
The flat
transition operator for driving a Markov chain over a space
of ensembles.
createSystemRandom :: IO GenIO
Seed a PRNG with data from the system's fast source of pseudo-random
numbers. All the caveats of withSystemRandom
apply here as well.
withSystemRandom :: PrimBase m => (Gen (PrimState m) -> m a) -> IO a
Seed a PRNG with data from the system's fast source of
pseudo-random numbers ("/dev/urandom
" on Unix-like systems or
RtlGenRandom
on Windows), then run the given action.
This is a somewhat expensive function, and is intended to be called
only occasionally (e.g. once per thread). You should use the Gen
it creates to generate many random numbers.