flat-mcmc-1.4.1: Painless general-purpose sampling.

Copyright (c) 2016 Jared Tobin MIT Jared Tobin unstable ghc None Haskell2010

Numeric.MCMC.Flat

Description

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

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.Unboxed (toList)
>>> :{
>>> let rosenbrock xs = negate (5  *(x1 - x0 ^ 2) ^ 2 + 0.05 * (1 - x0) ^ 2)
            where [x0, x1] = toList xs
>>> :}
>>> :{
>>> 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]
>>> ]
>>> :}
>>> withSystemRandom . asGenIO \$ mcmc 2 origin 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


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

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.

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.

create :: PrimMonad m => m (Gen (PrimState m)) #

Create a generator for variates using a fixed seed.

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.

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

Constrain the type of an action to run in the IO monad.

A type-specialized alias for Data.Vector.fromList.

Use this to create ensembles from lists of particles.

A type-specialized alias for Data.Vector.Unboxed.fromList

Use this to create particles from lists of doubles.