MCMC Metropolis–Hastings algorithm

Question

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][[2]] > 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

Answer 1

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.

Answer 2

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

ClearAll["Global`*"]

SeedRandom[123]
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, 
   ChartElementFunction -> "FadingRectangle",
   ChartStyle -> Orange], 
  Histogram[sim[[All, 1]], Axes -> {True, False}, ImageSize -> 300, 
   PlotLabel -> "b", PlotRange -> All, 
   ChartElementFunction -> "FadingRectangle",
   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}]