mcmc-0.5.0.0: Sample from a posterior using Markov chain Monte Carlo
Copyright(c) Dominik Schrempf 2021
LicenseGPL-3.0-or-later
Maintainerdominik.schrempf@gmail.com
Stabilityunstable
Portabilityportable
Safe HaskellNone
LanguageHaskell2010

Mcmc.Proposal

Description

Creation date: Wed May 20 13:42:53 2020.

Synopsis

Proposals and types

newtype PName Source #

Proposal name.

Constructors

PName 

Fields

Instances

Instances details
Eq PName Source # 
Instance details

Defined in Mcmc.Proposal

Methods

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

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

Ord PName Source # 
Instance details

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 Source # 
Instance details

Defined in Mcmc.Proposal

Methods

showsPrec :: Int -> PName -> ShowS #

show :: PName -> String #

showList :: [PName] -> ShowS #

Semigroup PName Source # 
Instance details

Defined in Mcmc.Proposal

Methods

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

sconcat :: NonEmpty PName -> PName #

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

Monoid PName Source # 
Instance details

Defined in Mcmc.Proposal

Methods

mempty :: PName #

mappend :: PName -> PName -> PName #

mconcat :: [PName] -> PName #

data PWeight Source #

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

Instances

Instances details
Eq PWeight Source # 
Instance details

Defined in Mcmc.Proposal

Methods

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

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

Ord PWeight Source # 
Instance details

Defined in Mcmc.Proposal

Show PWeight Source # 
Instance details

Defined in Mcmc.Proposal

pWeight :: Int -> PWeight Source #

Check if the weight is positive.

data PDimension Source #

Proposal dimension.

The number of affected, independent parameters.

The optimal acceptance rate of low dimensional proposals is higher than for high dimensional ones.

Optimal acceptance rates are still subject to controversies. As far as I know, research has focused on random walk proposal with a multivariate normal distribution of dimension d. In this case, the following acceptance rates are desired:

  • one dimension: 0.44 (numerical results);
  • five and more dimensions: 0.234 (numerical results);
  • infinite dimensions: 0.234 (theorem for specific target distributions).

See Handbook of Markov chain Monte Carlo, chapter 4.

Of course, many proposals may not be classical random walk proposals. For example, the beta proposal on a simplex (beta) samples one new variable of the simplex from a beta distribution while rescaling all other variables. What is the dimension of this proposal? I don't know, but I set the dimension to 2. The reason is that if the dimension of the simplex is 2, two variables are changed. If the dimension of the simplex is high, one variable is changed substantially, while all others are changed marginally.

Further, if a proposal changes a number of variables in the same way (and not independently like in a random walk proposal), I still set the dimension of the proposal to the number of variables changed.

Finally, I assume that proposals of unknown dimension have high dimension, and use the optimal acceptance rate 0.234.

data Proposal a Source #

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.

Predefined proposals are provided. To create custom proposals, one may use the convenience function createProposal.

Constructors

Proposal 

Fields

Instances

Instances details
Eq (Proposal a) Source # 
Instance details

Defined in Mcmc.Proposal

Methods

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

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

Ord (Proposal a) Source # 
Instance details

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 #

type KernelRatio = Log Double Source #

Ratio of the proposal kernels.

Part of the MHG acceptance ratio.

See also Jacobian.

NOTE: Actually the Jacobian should be part of the KernelRatio. However, it is more declarative to have them separate. It is a constant reminder: Is the Jacobian modifier different from 1.0?

type Jacobian = Log Double Source #

Absolute value of the determinant of the Jacobian matrix.

Part of the MHG acceptance ratio.

See also Jacobian.

type JacobianFunction a = a -> Jacobian Source #

Function calculating the Jacobian of a proposal.

(@~) :: Lens' b a -> Proposal a -> Proposal b infixl 7 Source #

Lift a proposal from one data type to another.

Assume the Jacobian is 1.0 (see also liftProposal and liftProposalWith).

For example:

scaleFirstEntryOfTuple = _1 @~ scale

liftProposal :: Lens' b a -> Proposal a -> Proposal b Source #

Lift a proposal from one data type to another.

Assume the Jacobian is 1.0 (see also (@~) and liftProposalWith).

liftProposalWith :: JacobianFunction b -> Lens' b a -> Proposal a -> Proposal b Source #

Lift a proposal from one data type to another.

A function to calculate the Jacobian has to be provided (see also liftProposal).

For further reference, please see the example Pair.

type ProposalSimple a = a -> GenIO -> IO (a, KernelRatio, Jacobian) Source #

Simple proposal without tuning information.

Instruction about randomly moving from the current state to a new state, given some source of randomness.

In order to calculate the Metropolis-Hastings-Green ratio, we need to know the ratio of the backward to forward kernels (the KernelRatio or the probability masses or probability densities) and the Jacobian.

For unbiased proposals, these values are 1.0 such that

proposalSimpleUnbiased x g = return (x', 1.0, 1.0)

For biased proposals, the kernel ratio is qYX / qXY, where qXY is the probability density to move from X to Y, and the absolute value of the determinant of the Jacobian matrix differs from 1.0.

data Tuner a Source #

Tune the acceptance rate of a Proposal; see tune, or autoTuneCycle.

data Tune Source #

Tune the proposal?

Constructors

Tune 
NoTune 

Instances

Instances details
Eq Tune Source # 
Instance details

Defined in Mcmc.Proposal

Methods

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

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

Show Tune Source # 
Instance details

Defined in Mcmc.Proposal

Methods

showsPrec :: Int -> Tune -> ShowS #

show :: Tune -> String #

showList :: [Tune] -> ShowS #

createProposal Source #

Arguments

:: PDescription

Description of the proposal type and parameters.

-> (TuningParameter -> ProposalSimple a)

Function creating a simple proposal for a given tuning parameter. The larger the tuning parameter, the larger the proposal and the lower the expected acceptance rate; and vice versa.

-> PDimension

Dimension.

-> PName

Name.

-> PWeight

Weight.

-> Tune

Activate tuning?

-> Proposal a 

Create a tuneable proposal.

type TuningParameter = Double Source #

Tuning parameter.

tuningParameterMin :: TuningParameter Source #

Minimal tuning parameter; 1e-12, subject to change.

>>> tuningParameterMin
1e-5

tuningParameterMax :: TuningParameter Source #

Maximal tuning parameter; 1e12, subject to change. >>> tuningParameterMax 1e3

tune :: (TuningParameter -> TuningParameter) -> Proposal a -> Maybe (Proposal a) Source #

Tune a Proposal.

The size of the proposal is proportional to the tuning parameter which has a positive lower bound of tuningParameterMin.

The tuning function maps the current tuning parameter to a new one.

Return Nothing if Proposal is not tuneable.

Cycles

data Order Source #

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

Instances details
Eq Order Source # 
Instance details

Defined in Mcmc.Proposal

Methods

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

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

Show Order Source # 
Instance details

Defined in Mcmc.Proposal

Methods

showsPrec :: Int -> Order -> ShowS #

show :: Order -> String #

showList :: [Order] -> ShowS #

Default Order Source # 
Instance details

Defined in Mcmc.Proposal

Methods

def :: Order #

data Cycle a Source #

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 Source #

Create a Cycle from a list of Proposals.

setOrder :: Order -> Cycle a -> Cycle a Source #

Set the order of Proposals in a Cycle.

prepareProposals :: Cycle a -> GenIO -> IO [Proposal a] Source #

Replicate Proposals according to their weights and possibly shuffle them.

autoTuneCycle :: Acceptance (Proposal a) -> Cycle a -> Cycle a Source #

Calculate acceptance rates and auto tune the Proposals in the Cycle. For now, a Proposal is enlarged when the acceptance rate is above 0.44, and shrunk otherwise. Do not change Proposals that are not tuneable.

Acceptance rates

data Acceptance k Source #

For each key k, store the number of accepted and rejected proposals.

Instances

Instances details
Eq k => Eq (Acceptance k) Source # 
Instance details

Defined in Mcmc.Proposal

Methods

(==) :: Acceptance k -> Acceptance k -> Bool #

(/=) :: Acceptance k -> Acceptance k -> Bool #

(Ord k, Read k) => Read (Acceptance k) Source # 
Instance details

Defined in Mcmc.Proposal

Show k => Show (Acceptance k) Source # 
Instance details

Defined in Mcmc.Proposal

ToJSONKey k => ToJSON (Acceptance k) Source # 
Instance details

Defined in Mcmc.Proposal

(Ord k, FromJSONKey k) => FromJSON (Acceptance k) Source # 
Instance details

Defined in Mcmc.Proposal

emptyA :: Ord k => [k] -> Acceptance k Source #

In the beginning there was the Word.

Initialize an empty storage of accepted/rejected values.

pushA :: Ord k => k -> Bool -> Acceptance k -> Acceptance k Source #

For key k, prepend an accepted (True) or rejected (False) proposal.

resetA :: Ord k => Acceptance k -> Acceptance k Source #

Reset acceptance storage.

transformKeysA :: (Ord k1, Ord k2) => [k1] -> [k2] -> Acceptance k1 -> Acceptance k2 Source #

Transform keys using the given lists. Keys not provided will not be present in the new Acceptance variable.

acceptanceRate :: Ord k => k -> Acceptance k -> Maybe (Int, Int, Double) Source #

Acceptance counts and rate for a specific proposal.

Return Nothing if no proposals have been accepted or rejected (division by zero).

acceptanceRates :: Acceptance k -> Map k (Maybe Double) Source #

Acceptance rates for all proposals.

Set rate to Nothing if no proposals have been accepted or rejected (division by zero).

Output

proposalHeader :: ByteString Source #

Header of proposal summaries.

proposalHLine :: ByteString Source #

Horizontal line of proposal summaries.

summarizeCycle :: Acceptance (Proposal a) -> Cycle a -> ByteString Source #

Summarize the Proposals in the Cycle. Also report acceptance rates.