Hi All, I am trying to implement Metropolis–Hastings algorithm to find parameters. I basically followed the pdf attachment. 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 https://mathematica.stackexchange.com/questions/152789/mcmc-metropolis-hastings-algorithm

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

Attachments: