Numeric.MCMC.Metropolis

This implementation uses spherical Gaussian proposals to implement a reliable and computationally inexpensive sampling routine. It can be used as a baseline from which to benchmark other algorithms for a given problem.

The mcmc function streams a trace to stdout to be processed elsewhere, while the metropolis transition can be used for more flexible purposes, such as working with samples in memory.

Documentation

mcmc :: (MonadIO m, PrimMonad m, Traversable f, Show (f Double)) => Int -> Double -> f Double -> (f Double -> Double) -> Gen (PrimState m) -> m ()

Trace n iterations of a Markov chain and stream them to stdout.

>>> let rosenbrock [x0, x1] = negate (5  *(x1 - x0 ^ 2) ^ 2 + 0.05 * (1 - x0) ^ 2)
>>> withSystemRandom . asGenIO $ mcmc 3 1 [0, 0] rosenbrock
0.5000462419822702,0.5693944056267897
0.5000462419822702,0.5693944056267897
-0.7525995304580824,1.2240725505283248

chain

Arguments

:: (PrimMonad m, Traversable f)

=> Int

Number of iterations

-> Double

Step standard deviation

-> f Double

Starting position

-> (f Double -> Double)

Log-density (to additive constant)

-> Gen (PrimState m)

PRNG

-> m [Chain (f Double) b]

Trace n iterations of a Markov chain and collect the results in a list.

>>> let rosenbrock [x0, x1] = negate (5  *(x1 - x0 ^ 2) ^ 2 + 0.05 * (1 - x0) ^ 2)
>>> results <- withSystemRandom . asGenIO $ chain 3 1 [0, 0] rosenbrock
>>> mapM_ print results
0.0,0.0
1.4754117657794871e-2,0.5033208261760778
3.8379699517007895e-3,0.24627131099479127

chain' :: (PrimMonad m, Traversable f) => Int -> Double -> f Double -> (f Double -> Double) -> Maybe (f Double -> b) -> Gen (PrimState m) -> m [Chain (f Double) b]

Return a list of Chain values potentially with tunable values computed from each position.

metropolis :: (Traversable f, PrimMonad m) => Double -> Maybe (f Double -> b) -> Transition m (Chain (f Double) b)

A generic Metropolis transition operator.

Re-exported

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

A generic transition operator.

Has access to randomness via the underlying Prob monad.

data Chain a b

The Chain type specifies the state of a Markov chain at any given iteration.

Constructors

Chain

Fields

• chainTarget :: Target a
• chainScore :: !Double
• chainPosition :: a
• chainTunables :: Maybe b

Instances

Show a => Show (Chain a b)

Defined in Data.Sampling.Types

Methods

showsPrec :: Int -> Chain a b -> ShowS

show :: Chain a b -> String

showList :: [Chain a b] -> ShowS

data Target a

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.

Constructors

Target

Fields

• lTarget :: a -> Double
• glTarget :: Maybe (a -> a)

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

Create a generator for variates using a fixed seed.

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.

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

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