29

These days I'm trying to conduct a model sensitivity test which is heavily based on the Markov Chain Monte Carlo simulation approach.

And I find this 'MCMC' package that can perform Markov Chain Monte Carlo simulations.

However, I found this package doesn't use much of the built-in stochastic process functions.

Also, due to my very limited knowledge of stochastic processes, I'd like to get better understanding of MCMC via a more Mathematica way.

As such, I'm wondering if such MCMC simulations can be performed with higher-level built-in stochastic functions (e.g. RandomFunction, HiddenMarkovProcess)?

Update on 20160330:

Now I get better understanding of how the MCMC procedure works.

But, one thing I'm still not clear is how to generate the transition matrix needed by the DsicreteMarkovProcess in Mathematica.

Can someone provide a clue for this?

Note: Cross-posted at http://community.wolfram.com/groups/-/m/t/830659

gwr
  • 13,452
  • 2
  • 47
  • 78
sunt05
  • 4,367
  • 23
  • 34
  • Thank you for noting that this is cross posted. Can you also note that on the community website as well? – Searke Mar 28 '16 at 15:03
  • Noted on the community. Thanks for the reminder, @Searke :) – sunt05 Mar 28 '16 at 15:30
  • 5
    The book "Bayesian Logical Data Analysis for the Physical Sciences", by Phil Gregory, comes with a MCMC program in Mathematica. Look under Resources at http://bit.ly/1tdaYhG . – Romke Bontekoe Mar 30 '16 at 09:50
  • 1
    The transition matrix is holds the probability of moving to another state from the current state. The row index is the current state and the column indices are the states that can be transitioned to. The values are the probabilities of transferring to a state. Therefore the sum of the rows should equal one and there should be no negative entries. (E.g. if row 2 = {.3, .1, .6} then when in state 2 there is .3 probability of moving to state 1, .1 of staying in 2 and .6 of moving to 3. – Edmund Mar 30 '16 at 13:13

2 Answers2

23

As of Version 11.3, there is an undocumented utility package, called 

Statistics`MCMC`

I came to know of it from the example notebook of this repository by @Sjoerd Smit. The main function to create the Markov Chain is 

Statistics`MCMC`BuildMarkovChain["method"][arguments]

Run

Statistics`MCMC`MCMCData[]

to get all available methods to run. To check the usage, put "Usage" after the name of the method. For example, to get the usage of the {"AdaptiveMetropolis", "Log"}, run:

Statistics`MCMC`MCMCData[{"AdaptiveMetropolis", "Log"}, "Usage"]  
Statistics`MCMC`MCMCData[{"AdaptiveMetropolis", "Log"}, "Example"]

Once the chain is created, use

Statistics`MCMC`MarkovChainIterate[chainName, numOfSamples]
Statistics`MCMC`MarkovChainIterate[chainName, {numOfSamples, numberOfIterationsBetweenEachSample}]

to iterate the chain. The output is the sample. To get the status of the chain, print it. The chain is a MarkovChainObject. Take its first part to get an Association. It contains all info of the chain at that point in the iteration.

Example (update)

(Following the example in the aforementioned example notebook)
Let's create a custom logPDF function, set an initial point, an initial covariance matrix, and the number of steps after which to calculate the covariance matrices. This lets us create a Markov chain.

logPDF=Compile[{{pt, _Real, 1}}, - pt. pt];
initPoint={0,0};
initCovariance=DiagonalMatrix[ConstantArray[1, Length[initPoint]]];
delay=20;
chain=Statistics`MCMC`BuildMarkovChain[{"AdaptiveMetropolis", "Log"}][initPoint,logPDF,{initCovariance,delay}, Real, Compiled -> True]

Out:=

Out:=

Now let's do a burn-in for 10000 steps, and then create a sample of 100000.

Statistics`MCMC`MarkovChainIterate[chain,10000];
sample=Statistics`MCMC`MarkovChainIterate[chain,100000];

Mean[sample]  
Covariance[sample]
Correlation[sample]
SmoothDensityHistogram[sample]

Out:=
{-0.0124793, -0.0125073}
{{0.526805, -0.00857408}, {-0.00857408, 0.496705}}
{{1., -0.0167615}, {-0.0167615, 1.}}
enter image description here

raja
  • 346
  • 2
  • 7
5

I'm currently implementing a MCMC in Mathematica. What I've done so far.

1.-Build your Model

2.-Import your data (X,Y,σ)

3.-Assign the initial parameters (P1=1 and P2=2) for example.

4.-Sample new proposed values for P1 and P2 (we call them pP1 and pP2). This is done by creating a Multinormal Distribution with means {P1,P2} and a squared covariance matrix. We'll use an identity matrix here {{1,0},{0,1}}. We used a normal distributions because they are symmetric

{pP1, pP2} = RandomVariate[MultinormalDistribution[{P1, P2}, {{1, 0}, {0, 1}}]]

5.- We evaluate the likelihood of our model p(D/M) with both sets of parameters. For this we compare the prediction of our models (M[Xi]) for each data point (Xi,Yi,σi):

$\frac{e^{\frac{-(\text{Yi}-M[\text{Xi}])}{2 \text{$\sigma \ $i}^2}}}{\sqrt{2 \pi } \text{$\sigma $i}}$

We multiply this for every data point and obtain the likelihood with the old p(D/M1) and new parameters p(D/M2)

6.- Then we calculate the Metroplis Ratio 

metropolisr = p(D/M1)/p(D/M2)

If the metropolisr ≥ 1 we accept the new parameters, if metropolisr < 1 we accept it with a probability equal to metropolisr. You can write this as:

If [metropolisr ≥ RandomReal[],{P1,P2} = {pP1,pP2}]

This is the basic approach I'm using. You just have to do the iterations next.

BPinto
  • 1,022
  • 6
  • 19