# MCMC Metropolis–Hastings algorithm

I am trying to implement Metropolis–Hastings algorithm to find parameters. I basically followed this PDF. It seems it works but gives me off parameters. Note: This is the Metropolis algorithm not Metropolis–Hastings algorithm. I know there is MCMC package is available but I want to understand method. Can someone help me out to modified it? Thanks in advance.

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

a0 = 2; (*parameter 1*)
b0 = 5;  (*parameter 2*)
y[x_] := a0 x + b0 (*model*)

data = y@Range@10;

prior = PDF[NormalDistribution[0, 1], a] PDF[NormalDistribution[0, 1], b]
likelihood = Likelihood[NormalDistribution[a, b], data]

posterior[a_, b_] = prior *likelihood

ratio[{a1_, b1_}, {a2_, b2_}] =  posterior[a2, b2]/posterior[a1, b1] // Simplify

proposal[\[Theta]_,  std_] := \[Theta] + (RandomReal[NormalDistribution[0, #]] & /@ std)

test[\[Phi]_] := If[\[Phi][] > 0, True, False]

metroPolisStep[\[Theta]_, std_] :=With[{\[Phi] = proposal[\[Theta], std], r = RandomReal[]},If[test[\[Phi]],  If[r <= ratio[\[Theta], \[Phi]], \[Phi], \[Theta]], \[Theta]]]

metroPolis[initialState_, std_, steps_] :=  NestList[metroPolisStep[#, std] &, initialState, steps]

numStep = 10000;

burnin = Ceiling[numStep* 10/100];

metro = metroPolis[{2, 1}, {0.5, 1}, numStep];

param1 = Drop[metro[[All, 1]], burnin];

param2 = Drop[metro[[All, 2]], burnin];

ListLinePlot[param1, PlotStyle -> PointSize[Tiny], AspectRatio -> 0.2,  ImageSize -> 400, PlotRange -> All]

ListLinePlot[param2, PlotStyle -> PointSize[Tiny], AspectRatio -> 0.2, ImageSize -> 400, PlotRange -> All]

{SmoothHistogram[param1, Axes -> {True, False}, ImageSize -> 300,    PlotLabel -> "a", PlotRange -> All],   SmoothHistogram[param2, Axes -> {True, False}, ImageSize -> 300,    PlotLabel -> "b", PlotRange -> All]} // GraphicsRow

• Since SE does not support attachments, you could have at least edited your question here to link to the PDF instead of just doing copy-pasta. – J. M. will be back soon Aug 1 '17 at 5:32
• Sorry, I did that. – OkkesDulgerci Aug 1 '17 at 15:11
• I'm confused by your text. You want to implement the Metropolis-Hastings algorithm but you state "Note: This is the Metropolis algorithm not Metropolis–Hastings algorithm." – JimB Aug 1 '17 at 15:20
• Main difference is Metropolis algorithm assumes dist. is symmetric and Metropolis-Hastings algorithm is more general case. – OkkesDulgerci Aug 1 '17 at 15:27
• Would this help? 12000.org/my_notes/hastings_metropolis/index.htm Note there is a link for download of a notebook. – bill s Aug 1 '17 at 17:55

The following does not answer the OP's question directly, in that it does not provide modifications of the code presented. The purpose of this "answer" is to provide a clear statement of the Metropolis-Hastings algorithm and its relation to the Metropolis algorithm in hopes that this would aid the OP in modifying the code him- or herself. My impression from reading the OP's question is that there may be a misunderstanding on the OP's part, but perhaps I am wrong and this "answer" should be removed.

Given data (vector) $y$ and unknown parameter (vector) $\theta$, the posterior distribution is given by $$p(\theta|y) = \frac{p(y|\theta)\,p(\theta)}{p(y)} ,$$ where $p(y|\theta)$ is the likelihood, $p(\theta)$ is the prior, and $p(y)$ is (often called) the marginal likelihood.

The Metropolis-Hastings algorithm can be summarized as follows. Let $q(\theta'|\theta)$ denote the density for the proposal distribution for $\theta'$ conditioned on $\theta$. Let $\theta^{(r)}$ denote the current state of the MCMC chain. Then \begin{equation} \theta^{(r+1)} = \begin{cases} \theta' & \mathcal{R} \ge u \\ \theta^{(r)} & \text{otherwise} \end{cases} , \end{equation} where $u \sim \textsf{Uniform}(0,1)$ and \begin{equation} \mathcal{R} = \frac{p(\theta'|y)}{p(\theta^{(r)}|y)} \times \frac{q(\theta^{(r)}|\theta')}{q(\theta'|\theta^{(r)})} . \end{equation} Note that $p(y)$ cancels out of the first factor on the right-hand side. Also note that if \begin{equation} q(\theta|\theta') = q(\theta'|\theta) \qquad\text{for all $(\theta,\theta')$} \end{equation} then the second factor disappears and we have the Metropolis algorithm.

• This doesn't seem to have anything to do with Mathematica. The OP wanted help with modifying the algorithm. A Mathematica implementation of what you wrote would be more in line with the question. – b3m2a1 Aug 14 '17 at 18:27
• @mef Thanks for your post. I have read about theory but I need help for implementation. – OkkesDulgerci Aug 16 '17 at 3:43

I tried this but I am still not sure if I am doing right.

ClearAll["Global*"]

SeedRandom
m0 = 3;
b0 = 5;
std = 5;
sampleSize = 31;
y[x_] := m0 x + b0

data = Thread[{Range[0, 30],
y@Range[0, 30] +
RandomVariate[NormalDistribution[0, std], sampleSize]}]

Show[Plot[y[x], {x, 0, 30}, Frame -> True, Axes -> False],
ListPlot[data, PlotStyle -> Red]]

model[x_] := m x + b

prior = PDF[NormalDistribution[3, 1], m] PDF[NormalDistribution[5, 1],
b]

Variance[data[[All, 2]]]

n = 31;
\[Sigma] = StandardDeviation[data[[All, 2]]]

likelihood =
Product[(1/Sqrt[2 \[Pi] (\[Sigma]^2) ]
Exp[-(1/(
2 \[Sigma]^2)) (data[[i, 2]] - model[data[[i, 1]]])^2]), {i, n}]

posterior[m_, b_] = prior likelihood

pdf = Compile[{{m, _Real}, {b, _Real}}, posterior[m, b],
Parallelization -> True];

mcmcfunComp =
Compile[{{s, _Real, 1}},
Module[{xm, ym, proposal, xp, yp, p2, p1, proposalSigma = 0.2},
{xm, ym} = s;
xp = RandomReal[NormalDistribution[xm, proposalSigma]];
yp = RandomReal[NormalDistribution[ym, proposalSigma]];
p2 = pdf[xp, yp];
p1 = pdf[xm, ym];
proposal = p2/p1;
If[RandomReal[] <= proposal, {xp, yp}, {xm, ym}]],
CompilationOptions -> {"InlineExternalDefinitions" -> True},
Parallelization -> True];

sim = Drop[NestList[mcmcfunComp[#] &, {2, 4}, 100000], 10000];

ListLinePlot[#, PlotRange -> All, PlotStyle -> Thin,
AspectRatio -> 0.3, Frame -> True,
ImageSize -> Medium] & /@ {sim[[All, 1]],
sim[[All, 2]]} // GraphicsRow

{Histogram[sim[[All, 1]], Axes -> {True, False}, ImageSize -> 300,
PlotLabel -> "a", PlotRange -> All,
ChartStyle -> Orange],
Histogram[sim[[All, 1]], Axes -> {True, False}, ImageSize -> 300,
PlotLabel -> "b", PlotRange -> All,
ChartStyle -> Orange]} // GraphicsRow

{Mean@sim[[All, 1]], Mean@sim[[All, 2]]}

{SmoothHistogram[sim[[All, 1]], Axes -> {True, False},
ImageSize -> 300, PlotLabel -> "m", PlotRange -> All],
SmoothHistogram[sim[[All, 2]], Axes -> {True, False},
ImageSize -> 300, PlotLabel -> "b",
PlotRange -> All]} // GraphicsRow

Labeled[SmoothDensityHistogram[Thread[{sim[[All, 1]], sim[[All, 2]]}],
"Oversmooth", "PDF", ColorFunction -> "ThermometerColors",
ImageSize -> 300, Mesh -> 3, MeshStyle -> Directive[White, Dashed],
GridLines -> Automatic,
PlotRange -> All], {Style["m ", 20, FontFamily -> "Times New Roman",
Bold], Style["b ", 20, FontFamily -> "Times New Roman",
Bold]}, {Bottom, Left}]
`