How to implement Markov Chain Monte Carlo with built-in functions?

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:=

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.}}

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.