Mcmc

Creation date: Tue May 5 18:01:15 2020.

For an introduction to Markov chain Monte Carlo (MCMC) samplers and update mechanisms using the Metropolis-Hastings-Green algorithm, please see Geyer, C. J., (2011), Introduction to Markov Chain Monte Carlo, In Handbook of Markov Chain Monte Carlo (pp. 45), CRC press.

This library focusses on classical Markov chain Monte Carlo algorithms such as the Metropolis-Hastings-Green [1] algorithm, or population methods involving parallel chains such as the Metropolic-coupled Markov chain Monte Carlo [2] algorithm. In particular, sequential Monte Carlo [3] algorithms following a moving posterior distribution are not provided.

An MCMC sampler can be run with mcmc, for example using the Metropolis-Hastings-Green algorithm mhg.

Usually, it is best to start with an example:

• Basic inference of the accuracy of an archer
• More involved examples

The import of this module alone should cover most use cases.

[1] Geyer, C. J. (2011), Introduction to markov chain monte carlo, In Handbook of Markov Chain Monte Carlo (pp. 45), CRC press.

[2] Geyer, C. J. (1991), Markov chain monte carlo maximum likelihood, Computing Science and Statistics, Proceedings of the 23rd Symposium on the Interface.

[3] Sequential monte carlo methods in practice (2001), Editors: Arnaud Doucet, Nando de Freitas, and Neil Gordon, Springer New York.

Proposals

A Proposal is an instruction about how to advance a given Markov chain so that it possibly reaches a new state. That is, Proposals specify how the chain traverses the state space. As far as this MCMC library is concerned, Proposals are considered to be /elementary updates/ in that they cannot be decomposed into smaller updates.

Proposals can be combined to form composite updates, a technique often referred to as composition. On the other hand, mixing (used in the sense of mixture models) is the random choice of a Proposal (or a composition of Proposals) from a given set.

The composition and mixture of Proposals allows specification of nearly all MCMC algorithms involving a single chain (i.e., population methods such as particle filters are excluded). In particular, Gibbs samplers of all sorts can be specified using this procedure. For reference, please see the short encyclopedia of MCMC methods.

This library enables composition and mixture of Proposals via the Cycle data type. Essentially, a Cycle is a set of Proposals. The chain advances after the completion of each Cycle, which is called an iteration, and the iteration counter is increased by one.

The Proposals in a Cycle can be executed in the given order or in a random sequence which allows, for example, specification of a fixed scan Gibbs sampler, or a random sequence scan Gibbs sampler, respectively. See Order.

Note that it is of utter importance that the given Cycle enables traversal of the complete state space. Otherwise, the Markov chain will not converge to the correct stationary posterior distribution.

Proposals are named according to what they do, i.e., how they change the state of a Markov chain, and not according to the intrinsically used probability distributions. For example, slideSymmetric is a sliding proposal. Under the hood, it uses the normal distribution with mean zero and given variance. The sampled variate is added to the current value of the variable (hence, the name slide). The same nomenclature is used by RevBayes [4]. The probability distributions and intrinsic properties of a specific proposal are specified in detail in the documentation.

The other method, which is used intrinsically, is more systematic, but also a little bit more complicated: we separate between the proposal distribution and how the state is affected. And here, I am referring to the operator (addition, multiplication, any other binary operator). For example, the sliding proposal with mean m, standard deviation s, and tuning parameter t is implemented as

slideSimple :: Double -> Double -> Double -> ProposalSimple Double
slideSimple m s t =
  genericContinuous (normalDistr m (s * t)) (+) (Just negate) Nothing

This specification is more involved. Especially since we need to know the probability of jumping back, and so we need to know the inverse operator negate. However, it also allows specification of new proposals with great ease.

[4] Höhna, S., Landis, M. J., Heath, T. A., Boussau, B., Lartillot, N., Moore, B. R., Huelsenbeck, J. P., …, Revbayes: bayesian phylogenetic inference using graphical models and an interactive model-specification language, Systematic Biology, 65(4), 726–736 (2016). http://dx.doi.org/10.1093/sysbio/syw021

newtype PName

Proposal name.

Constructors

PName

Fields

• fromPName :: String

Instances

Eq PName

Defined in Mcmc.Proposal

Methods

(==) :: PName -> PName -> Bool

(/=) :: PName -> PName -> Bool

Ord PName

Defined in Mcmc.Proposal

Methods

compare :: PName -> PName -> Ordering

(<) :: PName -> PName -> Bool

(<=) :: PName -> PName -> Bool

(>) :: PName -> PName -> Bool

(>=) :: PName -> PName -> Bool

max :: PName -> PName -> PName

min :: PName -> PName -> PName

Show PName

Defined in Mcmc.Proposal

Methods

showsPrec :: Int -> PName -> ShowS

show :: PName -> String

showList :: [PName] -> ShowS

Semigroup PName

Defined in Mcmc.Proposal

Methods

(<>) :: PName -> PName -> PName

sconcat :: NonEmpty PName -> PName

stimes :: Integral b => b -> PName -> PName

Monoid PName

Defined in Mcmc.Proposal

Methods

mempty :: PName

mappend :: PName -> PName -> PName

mconcat :: [PName] -> PName

newtype PWeight

The weight determines how often a Proposal is executed per iteration of the Markov chain.

Constructors

PWeight

Fields

• fromPWeight :: Int

Instances

Eq PWeight

Defined in Mcmc.Proposal

Methods

(==) :: PWeight -> PWeight -> Bool

(/=) :: PWeight -> PWeight -> Bool

Ord PWeight

Defined in Mcmc.Proposal

Methods

compare :: PWeight -> PWeight -> Ordering

(<) :: PWeight -> PWeight -> Bool

(<=) :: PWeight -> PWeight -> Bool

(>) :: PWeight -> PWeight -> Bool

(>=) :: PWeight -> PWeight -> Bool

max :: PWeight -> PWeight -> PWeight

min :: PWeight -> PWeight -> PWeight

Show PWeight

Defined in Mcmc.Proposal

Methods

showsPrec :: Int -> PWeight -> ShowS

show :: PWeight -> String

showList :: [PWeight] -> ShowS

data Proposal a

A Proposal is an instruction about how the Markov chain will traverse the state space a. Essentially, it is a probability mass or probability density conditioned on the current state (i.e., a Markov kernel).

A Proposal may be tuneable in that it contains information about how to enlarge or shrink the step size to tune the acceptance rate.

Instances

Eq (Proposal a)

Defined in Mcmc.Proposal

Methods

(==) :: Proposal a -> Proposal a -> Bool

(/=) :: Proposal a -> Proposal a -> Bool

Ord (Proposal a)

Defined in Mcmc.Proposal

Methods

compare :: Proposal a -> Proposal a -> Ordering

(<) :: Proposal a -> Proposal a -> Bool

(<=) :: Proposal a -> Proposal a -> Bool

(>) :: Proposal a -> Proposal a -> Bool

(>=) :: Proposal a -> Proposal a -> Bool

max :: Proposal a -> Proposal a -> Proposal a

min :: Proposal a -> Proposal a -> Proposal a

(@~) :: Lens' b a -> Proposal a -> Proposal b

Convert a proposal from one data type to another using a lens.

For example:

scaleFirstEntryOfTuple = _1 @~ scale

data Tune

Tune the proposal?

Constructors

Tune

NoTune

Instances

Eq Tune

Defined in Mcmc.Proposal

Methods

(==) :: Tune -> Tune -> Bool

(/=) :: Tune -> Tune -> Bool

Show Tune

Defined in Mcmc.Proposal

Methods

showsPrec :: Int -> Tune -> ShowS

show :: Tune -> String

showList :: [Tune] -> ShowS

scale

Arguments

:: Double

Shape.

-> Double

Scale.

-> PName

-> PWeight

-> Tune

-> Proposal Double

Multiplicative proposal with gamma distributed kernel.

scaleUnbiased

Arguments

:: Double

Shape.

-> PName

-> PWeight

-> Tune

-> Proposal Double

Multiplicative proposal with gamma distributed kernel.

The scale of the gamma distribution is set to (shape)^{-1}, so that the mean of the gamma distribution is 1.0.

scaleContrarily

Arguments

:: Double

Shape.

-> Double

Scale.

-> PName

-> PWeight

-> Tune

-> Proposal (Double, Double)

Multiplicative proposal with gamma distributed kernel.

The two values are scaled contrarily so that their product stays constant. Contrary proposals are useful when parameters are confounded.

scaleBactrian

Arguments

:: Double

Spike parameter.

-> Double

Standard deviation.

-> PName

Name.

-> PWeight

Weight.

-> Tune

Enable tuning.

-> Proposal Double

Multiplicative proposal with kernel similar to the silhouette of a Bactrian camel. See slideBactrian.

slide

Arguments

:: Double

Mean.

-> Double

Standard deviation.

-> PName

Name.

-> PWeight

Weight.

-> Tune

Enable tuning.

-> Proposal Double

Additive proposal with normally distributed kernel.

slideSymmetric

Arguments

:: Double

Standard deviation.

-> PName

Name.

-> PWeight

Weight.

-> Tune

Enable tuning.

-> Proposal Double

Additive proposal with normally distributed kernel with mean zero. This proposal is very fast, because the Metropolis-Hastings-Green ratio does not include calculation of the forwards and backwards kernels.

slideUniformSymmetric

Arguments

:: Double

Delta.

-> PName

Name.

-> PWeight

Weight.

-> Tune

Enable tuning.

-> Proposal Double

Additive proposal with uniformly distributed kernel with mean zero. This proposal is very fast, because the Metropolis-Hastings-Green ratio does not include calculation of the forwards and backwards kernels.

slideContrarily

Arguments

:: Double

Mean.

-> Double

Standard deviation.

-> PName

Name.

-> PWeight

Weight.

-> Tune

Enable tuning.

-> Proposal (Double, Double)

Additive proposal with normally distributed kernel.

The two values are slid contrarily so that their sum stays constant. Contrary proposals are useful when parameters are confounded.

slideBactrian

Arguments

:: Double

Spike parameter \(m\).

-> Double

Standard deviation \(s\).

-> PName

Name.

-> PWeight

Weight.

-> Tune

Enable tuning.

-> Proposal Double

Additive symmetric proposal with kernel similar to the silhouette of a Bactrian camel.

The Bactrian kernel is a mixture of two symmetrically arranged normal distributions. The spike parameter \(m \in (0, 1)\) loosely determines the standard deviations of the individual humps while the second parameter \(s > 0\) refers to the standard deviation of the complete Bactrian kernel.

See https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3845170/.

module Mcmc.Proposal.Simplex

data Cycle a

In brief, a Cycle is a list of proposals.

The state of the Markov chain will be logged only after all Proposals in the Cycle have been completed, and the iteration counter will be increased by one. The order in which the Proposals are executed is specified by Order. The default is RandomO.

No proposals with the same name and description are allowed in a Cycle, so that they can be uniquely identified.

cycleFromList :: [Proposal a] -> Cycle a

Create a Cycle from a list of Proposals.

data Order

Define the order in which Proposals are executed in a Cycle. The total number of Proposals per Cycle may differ between Orders (e.g., compare RandomO and RandomReversibleO).

Constructors

RandomO

Shuffle the Proposals in the Cycle. The Proposals are replicated according to their weights and executed in random order. If a Proposal has weight w, it is executed exactly w times per iteration.

SequentialO

The Proposals are executed sequentially, in the order they appear in the Cycle. Proposals with weight w>1 are repeated immediately w times (and not appended to the end of the list).

RandomReversibleO

Similar to RandomO. However, a reversed copy of the list of shuffled Proposals is appended such that the resulting Markov chain is reversible. Note: the total number of Proposals executed per cycle is twice the number of RandomO.

SequentialReversibleO

Similar to SequentialO. However, a reversed copy of the list of sequentially ordered Proposals is appended such that the resulting Markov chain is reversible.

Instances

Eq Order

Defined in Mcmc.Proposal

Methods

(==) :: Order -> Order -> Bool

(/=) :: Order -> Order -> Bool

Show Order

Defined in Mcmc.Proposal

Methods

showsPrec :: Int -> Order -> ShowS

show :: Order -> String

showList :: [Order] -> ShowS

Default Order

Defined in Mcmc.Proposal

Methods

def :: Order

setOrder :: Order -> Cycle a -> Cycle a

Set the order of Proposals in a Cycle.

Settings

module Mcmc.Settings

Monitors

A Monitor describes which part of the Markov chain should be logged and where. Monitor files can directly be loaded into established MCMC analysis tools working with tab separated tables (for example, Tracer).

There are three different Monitor types:

MonitorStdOut

Log to standard output.

MonitorFile

Log to a file.

MonitorBatch

Log summary statistics such as the mean of the last states to a file.

data Monitor a

A Monitor observing the chain.

A Monitor describes which part of the Markov chain should be logged and where. Further, they allow output of summary statistics per iteration in a flexible way.

Constructors

Monitor (MonitorStdOut a) [MonitorFile a] [MonitorBatch a]

data MonitorStdOut a

Monitor to standard output; constructed with monitorStdOut.

monitorStdOut

Arguments

:: [MonitorParameter a]

Instructions about which parameters to log.

-> Int

Logging period.

-> MonitorStdOut a

Monitor to standard output.

data MonitorFile a

Monitor to a file; constructed with monitorFile.

monitorFile

Arguments

:: String

Name; used as part of the file name.

-> [MonitorParameter a]

Instructions about which parameters to log.

-> Int

Logging period.

-> MonitorFile a

Monitor parameters to a file.

data MonitorBatch a

Batch monitor to a file.

Calculate summary statistics over the last given number of iterations (batch size). Construct with monitorBatch.

monitorBatch

Arguments

:: String

Name; used as part of the file name.

-> [MonitorParameterBatch a]

Instructions about how to calculate the summary statistics.

-> Int

Batch size.

-> MonitorBatch a

Batch monitor parameters to a file, see MonitorBatch.

module Mcmc.Monitor.Parameter

module Mcmc.Monitor.ParameterBatch

Prior distributions

Convenience functions for computing priors.

module Mcmc.Prior

MCMC samplers

mcmc :: Algorithm a => Settings -> a -> IO a

Run an MCMC algorithm with given settings.

mcmcContinue :: Algorithm a => Int -> Settings -> a -> IO a

Continue an MCMC algorithm for the given number of iterations.

Currently, it is only possible to continue MCMC algorithms that have completed successfully. This restriction is necessary, because for parallel chains, it is hardly possible to ensure all chains are synchronized when the process is killed.

See:

• mhgLoad
• mc3Load

See also settingsLoad, mhgLoad, and mc3Load.

Algorithms

module Mcmc.Algorithm.Metropolis

module Mcmc.Algorithm.MC3

Useful type synonyms

type PriorFunction a = a -> Log Double

Prior function.

noPrior :: PriorFunction a

Flat prior function. Useful for testing and debugging.

type LikelihoodFunction a = a -> Log Double

Likelihood function.

noLikelihood :: LikelihoodFunction a

Flat likelihood function. Useful for testing and debugging.