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.

## Synopsis

- mcmc :: (MonadIO m, PrimMonad m) => Int -> Ensemble -> (Particle -> Double) -> Gen (PrimState m) -> m ()
- flat :: PrimMonad m => Transition m Chain
- type Particle = Vector Double
- type Ensemble = Vector Particle
- data Chain
- type Transition (m :: Type -> Type) a = StateT a (Prob m) ()
- data Target a = Target {}
- 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
- ensemble :: [Vector Double] -> Vector (Vector Double)
- particle :: [Double] -> Vector Double

# Documentation

mcmc :: (MonadIO m, PrimMonad m) => Int -> Ensemble -> (Particle -> Double) -> Gen (PrimState m) -> m () 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.Unboxed (toList)`

`>>>`

`:{`

`>>>`

where [x0, x1] = toList xs`let rosenbrock xs = negate (5 *(x1 - x0 ^ 2) ^ 2 + 0.05 * (1 - x0) ^ 2)`

`>>>`

`:}`

`>>>`

`:{`

`>>>`

`let origin = ensemble [`

`>>>`

`particle [negate 1.0, negate 1.0]`

`>>>`

`, particle [negate 1.0, 1.0]`

`>>>`

`, particle [1.0, negate 1.0]`

`>>>`

`, particle [1.0, 1.0]`

`>>>`

`]`

`>>>`

`:}`

`>>>`

-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`withSystemRandom . asGenIO $ mcmc 2 origin rosenbrock`

flat :: PrimMonad m => Transition m Chain Source #

The `flat`

transition operator for driving a Markov chain over a space
of ensembles.

type Particle = Vector Double Source #

A particle is an n-dimensional point in Euclidean space.

You can create a particle by using the `particle`

helper function, or just
use Data.Vector.Unboxed.fromList.

type Ensemble = Vector Particle Source #

An ensemble is a collection of particles.

The Markov chain we're interested in will run over the space of ensembles, so you'll want to build an ensemble out of a reasonable number of particles to kick off the chain.

You can create an ensemble by using the `ensemble`

helper function, or just
use Data.Vector.fromList.

type Transition (m :: Type -> Type) a = StateT a (Prob m) () #

A generic transition operator.

Has access to randomness via the underlying `Prob`

monad.

A `Target`

consists of a function from parameter space to the reals, as
well as possibly a gradient.

Most implementations assume a *log*-target, so records are named
accordingly.

createSystemRandom :: IO GenIO #

Seed a PRNG with data from the system's fast source of pseudo-random numbers.

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, then run the given action.

This function is unsafe and for example allows STRefs or any other mutable data structure to escape scope:

`>>>`

`ref <- withSystemRandom $ \_ -> newSTRef 1`

`>>>`

2`withSystemRandom $ \_ -> modifySTRef ref succ >> readSTRef ref`

`>>>`

3`withSystemRandom $ \_ -> modifySTRef ref succ >> readSTRef ref`

asGenIO :: (GenIO -> IO a) -> GenIO -> IO a #

Constrain the type of an action to run in the `IO`

monad.