declarative-0.5.1: DIY Markov Chains.

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

Numeric.MCMC

Contents

Description

This module presents a simple combinator language for Markov transition operators that are useful in MCMC.

Any transition operators sharing the same stationary distribution and obeying the Markov and reversibility properties can be combined in a couple of ways, such that the resulting operator preserves the stationary distribution and desirable properties amenable for MCMC.

We can deterministically concatenate operators end-to-end, or sample from a collection of them according to some probability distribution. See Geyer, 2005 for details.

The result is a simple grammar for building composite, property-preserving transition operators from existing ones:

transition ::= primitive transition
| concatT transition transition
| sampleT transition transition


In addition to the above, this module provides a number of combinators for building composite transition operators. It re-exports a number of production-quality transition operators from the mighty-metropolis, speedy-slice, and hasty-hamiltonian libraries.

Markov chains can then be run over arbitrary Targets using whatever transition operator is desired.

import Numeric.MCMC
import Data.Sampling.Types

target :: [Double] -> Double
target [x0, x1] = negate (5  *(x1 - x0 ^ 2) ^ 2 + 0.05 * (1 - x0) ^ 2)

rosenbrock :: Target [Double]
rosenbrock = Target target Nothing

transition :: Transition IO (Chain [Double] b)
transition =
concatT
(sampleT (metropolis 0.5) (metropolis 1.0))
(sampleT (slice 2.0) (slice 3.0))

main :: IO ()
main = withSystemRandom . asGenIO $mcmc 10000 [0, 0] transition rosenbrock See the attached test suite for other examples. Synopsis # Documentation concatT :: Monad m => Transition m a -> Transition m a -> Transition m a Source # Deterministically concat transition operators together. concatAllT :: Monad m => [Transition m a] -> Transition m a Source # Deterministically concat a list of transition operators together. sampleT :: PrimMonad m => Transition m a -> Transition m a -> Transition m a Source # Probabilistically concat transition operators together. sampleAllT :: PrimMonad m => [Transition m a] -> Transition m a Source # Probabilistically concat transition operators together via a uniform distribution. bernoulliT :: PrimMonad m => Double -> Transition m a -> Transition m a -> Transition m a Source # Probabilistically concat transition operators together using a Bernoulli distribution with the supplied success probability. This is just a generalization of sampleT. frequency :: PrimMonad m => [(Int, Transition m a)] -> Transition m a Source # Probabilistically concat transition operators together using the supplied frequency distribution. This function is more-or-less an exact copy of frequency, except here applied to transition operators. anneal :: (Monad m, Functor f) => Double -> Transition m (Chain (f Double) b) -> Transition m (Chain (f Double) b) Source # An annealing transformer. When executed, the supplied transition operator will execute over the parameter space annealed to the supplied inverse temperature. let annealedTransition = anneal 0.30 (slice 0.5) mcmc :: (MonadIO m, PrimMonad m, Show (t a)) => Int -> t a -> Transition m (Chain (t a) b) -> Target (t a) -> Gen (PrimState m) -> m () Source # Trace n iterations of a Markov chain and stream them to stdout. >>> withSystemRandom . asGenIO$ mcmc 3 [0, 0] (metropolis 0.5) rosenbrock
-0.48939312153007863,0.13290702689491818
1.4541485365128892e-2,-0.4859905564050404
0.22487398491619448,-0.29769783186855125


chain :: (MonadIO m, PrimMonad m) => Int -> t a -> Transition m (Chain (t a) b) -> Target (t a) -> Gen (PrimState m) -> m [Chain (t a) b] Source #

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

>>> results <- withSystemRandom . asGenIO \$ chain 3 [0, 0] (metropolis 0.5) rosenbrock


# Re-exported

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

A generic Metropolis transition operator.

hamiltonian :: (Num (IxValue (t Double)), Traversable t, FunctorWithIndex (Index (t Double)) t, Ixed (t Double), PrimMonad m, (~) * (IxValue (t Double)) Double) => Double -> Int -> Transition m (Chain (t Double) b) #

A Hamiltonian transition operator.

slice :: (PrimMonad m, FoldableWithIndex (Index (t a)) t, Ixed (t a), Num (IxValue (t a)), Variate (IxValue (t a))) => IxValue (t a) -> Transition m (Chain (t a) b) #

A slice sampling transition operator.

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.

class Monad m => PrimMonad m #

Class of monads which can perform primitive state-transformer actions

Minimal complete definition

primitive

Associated Types

type PrimState (m :: * -> *) :: * #

State token type

Instances

 Associated Typestype PrimState (IO :: * -> *) :: * # Methodsprimitive :: (State# (PrimState IO) -> (#VoidRep, PtrRepLifted, State# (PrimState IO), a#)) -> IO a # PrimMonad (ST s) Associated Typestype PrimState (ST s :: * -> *) :: * # Methodsprimitive :: (State# (PrimState (ST s)) -> (#VoidRep, PtrRepLifted, State# (PrimState (ST s)), a#)) -> ST s a # PrimMonad m => PrimMonad (Prob m) Associated Typestype PrimState (Prob m :: * -> *) :: * # Methodsprimitive :: (State# (PrimState (Prob m)) -> (#VoidRep, PtrRepLifted, State# (PrimState (Prob m)), a#)) -> Prob m a # PrimMonad m => PrimMonad (ListT m) Associated Typestype PrimState (ListT m :: * -> *) :: * # Methodsprimitive :: (State# (PrimState (ListT m)) -> (#VoidRep, PtrRepLifted, State# (PrimState (ListT m)), a#)) -> ListT m a # PrimMonad m => PrimMonad (MaybeT m) Associated Typestype PrimState (MaybeT m :: * -> *) :: * # Methodsprimitive :: (State# (PrimState (MaybeT m)) -> (#VoidRep, PtrRepLifted, State# (PrimState (MaybeT m)), a#)) -> MaybeT m a # PrimMonad m => PrimMonad (IdentityT * m) Associated Typestype PrimState (IdentityT * m :: * -> *) :: * # Methodsprimitive :: (State# (PrimState (IdentityT * m)) -> (#VoidRep, PtrRepLifted, State# (PrimState (IdentityT * m)), a#)) -> IdentityT * m a # (Error e, PrimMonad m) => PrimMonad (ErrorT e m) Associated Typestype PrimState (ErrorT e m :: * -> *) :: * # Methodsprimitive :: (State# (PrimState (ErrorT e m)) -> (#VoidRep, PtrRepLifted, State# (PrimState (ErrorT e m)), a#)) -> ErrorT e m a # PrimMonad m => PrimMonad (ExceptT e m) Associated Typestype PrimState (ExceptT e m :: * -> *) :: * # Methodsprimitive :: (State# (PrimState (ExceptT e m)) -> (#VoidRep, PtrRepLifted, State# (PrimState (ExceptT e m)), a#)) -> ExceptT e m a # PrimMonad m => PrimMonad (StateT s m) Associated Typestype PrimState (StateT s m :: * -> *) :: * # Methodsprimitive :: (State# (PrimState (StateT s m)) -> (#VoidRep, PtrRepLifted, State# (PrimState (StateT s m)), a#)) -> StateT s m a # PrimMonad m => PrimMonad (StateT s m) Associated Typestype PrimState (StateT s m :: * -> *) :: * # Methodsprimitive :: (State# (PrimState (StateT s m)) -> (#VoidRep, PtrRepLifted, State# (PrimState (StateT s m)), a#)) -> StateT s m a # (Monoid w, PrimMonad m) => PrimMonad (WriterT w m) Associated Typestype PrimState (WriterT w m :: * -> *) :: * # Methodsprimitive :: (State# (PrimState (WriterT w m)) -> (#VoidRep, PtrRepLifted, State# (PrimState (WriterT w m)), a#)) -> WriterT w m a # (Monoid w, PrimMonad m) => PrimMonad (WriterT w m) Associated Typestype PrimState (WriterT w m :: * -> *) :: * # Methodsprimitive :: (State# (PrimState (WriterT w m)) -> (#VoidRep, PtrRepLifted, State# (PrimState (WriterT w m)), a#)) -> WriterT w m a # PrimMonad m => PrimMonad (ReaderT * r m) Associated Typestype PrimState (ReaderT * r m :: * -> *) :: * # Methodsprimitive :: (State# (PrimState (ReaderT * r m)) -> (#VoidRep, PtrRepLifted, State# (PrimState (ReaderT * r m)), a#)) -> ReaderT * r m a # (Monoid w, PrimMonad m) => PrimMonad (RWST r w s m) Associated Typestype PrimState (RWST r w s m :: * -> *) :: * # Methodsprimitive :: (State# (PrimState (RWST r w s m)) -> (#VoidRep, PtrRepLifted, State# (PrimState (RWST r w s m)), a#)) -> RWST r w s m a # (Monoid w, PrimMonad m) => PrimMonad (RWST r w s m) Associated Typestype PrimState (RWST r w s m :: * -> *) :: * # Methodsprimitive :: (State# (PrimState (RWST r w s m)) -> (#VoidRep, PtrRepLifted, State# (PrimState (RWST r w s m)), a#)) -> RWST r w s m a #

type family PrimState (m :: * -> *) :: * #

State token type

Instances

 type PrimState IO type PrimState IO = RealWorld type PrimState (ST s) type PrimState (ST s) = s type PrimState (Prob m) type PrimState (Prob m) = PrimState m type PrimState (ListT m) type PrimState (ListT m) = PrimState m type PrimState (MaybeT m) type PrimState (MaybeT m) = PrimState m type PrimState (IdentityT * m) type PrimState (IdentityT * m) = PrimState m type PrimState (ErrorT e m) type PrimState (ErrorT e m) = PrimState m type PrimState (ExceptT e m) type PrimState (ExceptT e m) = PrimState m type PrimState (StateT s m) type PrimState (StateT s m) = PrimState m type PrimState (StateT s m) type PrimState (StateT s m) = PrimState m type PrimState (WriterT w m) type PrimState (WriterT w m) = PrimState m type PrimState (WriterT w m) type PrimState (WriterT w m) = PrimState m type PrimState (ReaderT * r m) type PrimState (ReaderT * r m) = PrimState m type PrimState (RWST r w s m) type PrimState (RWST r w s m) = PrimState m type PrimState (RWST r w s m) type PrimState (RWST r w s m) = PrimState m

data RealWorld :: * #

RealWorld is deeply magical. It is primitive, but it is not unlifted (hence ptrArg). We never manipulate values of type RealWorld; it's only used in the type system, to parameterise State#.