# Metropolis-Hastings Algorithm Problem

I've been trying to implement the Metropolis-Hastings algorithm for a while but there seems to be something weird going on.

This algorithm does not need as much of the statistics to understand, and has the step size h.

I am working with the interval $(-\infty,\infty)$, if the interval is $[a,b]$ the module will define the function to be zero outside $[a,b]$ using the UnitStep command.

Here's the Metropolis algorithm used:

1. Choose $x_{0}$ as starting point.
2. Choose a step size $h$.
3. Generate $x_{trial}$=RandomVariate[UniformDistribution[{$x_{0}$-$\dfrac{h}{2}$,$x_{0}$+$\dfrac{h}{2}$}],1][1].
4. r = w($x_{trial}$)/w($x_{0}$).
5. If r$\geq$1 then $x_{1}$=$x_{trial}$.
6. If r$\leq$1 then the trial is accepted with a probability r. Choose a random number $\eta \in [0,1]$, if $\eta<r$ then $x_{1}$=$x_{trial}$, else $x_{1}$=$x_{0}$.

7.Repeat.

The function I'm using is $\dfrac{0.8}{1 - e^{-0.8 \pi}}e^{-0.8 x}$. It's the $f$ in the modules below. The points are sampled from a uniform distribution defined in the interval $[0, \pi]$; in the module $a=0$ and $b=\pi$. There is literally one difference between the two codes and that's step 4. In the Metropolis-Hastings algorithm you have the extra part added in the second code block but in the Metropolis there isn't such a thing. The only reason why the Metropolis works for the function is because I have added a step function to make areas outside the interval of $[0,\pi]$ to be zero.

Now, for the weirdness. If I simply use the metropolis algorithm, I get the following histogram which looks great.

This histogram looks fine. Here's the code for this histogram (you can copy it below).

However, when I run the Metropolis-Hastings algorithm. I get the following histogram.

This histogram has slight overestimations and underestimations. If I run the code with T=100000. I get the following histogram where the difference is quite significant. This does not happen with the first code block.

Here's the code for this histogram (you can copy it below).

I don't get why they are so different - why doesn't the second histogram match the graph well? I did a lot of research and couldn't figure it out. Maybe it's something with my code. Any help would be appreciated.

Here's the code you can copy. Code 1:

MCHastings[f_, a_, b_, h_, X0_, T_] :=
Module[{w, X, data, xtrial, r, i, \[Eta], p1, p2},

w = (1 - UnitStep[a - x] - UnitStep[-b + x])*f;

p1 = Plot[w, {x, Max[a - 1, -10], Min[b + 1, 10]},
PlotStyle -> {Thick}, PlotRange -> All];

X = X0;(*Step 1*)
data = {X};
Do[
xtrial =
RandomVariate[UniformDistribution[{X - h/2, X + h/2}], 1][[
1]];(*Step 3*)
r = (w /. x -> xtrial)/(w /. x -> X)  ;(*Step 4*)

If[r >= 1,
X = xtrial, {\[Eta] = RandomReal[{0, 1}],
If[\[Eta] < r, X = xtrial, X = X]}];(*Steps 5 and 6*)
AppendTo[data, X];
, {i, 2, T}];(*Step 7*)

p2 = Histogram[data, 50, "PDF"];

Return[Show[p2, p1, BaseStyle -> Medium, AxesOrigin -> {0, 0}]];
]


Code 2:

MCHastings[f_, a_, b_, h_, X0_, T_] :=

Module[{w, X, data, xtrial, r, i, \[Eta], p1, p2},

w = (1 - UnitStep[a - x] - UnitStep[-b + x])*f;

p1 = Plot[w, {x, Max[a - 1, -10], Min[b + 1, 10]},
PlotStyle -> {Thick}, PlotRange -> All];

X = X0;(*Step 1*)
data = {X};
Do[
xtrial =
RandomVariate[UniformDistribution[{X - h/2, X + h/2}], 1][[
1]];(*Step 3*)
r = (w /. x -> xtrial)/(w /. x -> X) *
RandomVariate[UniformDistribution[{X - h/2, X + h/2}], 1][[1]]/
RandomVariate[UniformDistribution[{xtrial - h/2, xtrial + h/2}],
1][[1]];  (*Step 4*)

If[r >= 1,
X = xtrial, {\[Eta] = RandomReal[{0, 1}],
If[\[Eta] < r, X = xtrial, X = X]}];(*Steps 5 and 6*)
AppendTo[data, X];
, {i, 2, T}];(*Step 7*)

p2 = Histogram[data, 50, "PDF"];

Return[Show[p2, p1, BaseStyle -> Medium, AxesOrigin -> {0, 0}]];
]

• Why are the histograms different? It makes it harder to see the difference. Also, you could write the code in a way to avoid having the same instructions multiple times (it would also make it easier to spot the differences in the two functions). Feb 19, 2018 at 3:27
• Please explain the algorithms in plain English (in addition to the code). I.e. what problem you're applying it to, what space you're sampling from, what's your proposal distribution, etc. Don't make people figure this out from (potentially buggy) code when it could be stated in a few sentences. Feb 19, 2018 at 9:41
• You could try reproducing the algorithm in this wikipedia section. Feb 20, 2018 at 20:39
• @Szabolcs I've added comments above which should appear after peer review. Please let me know if it helps. Feb 21, 2018 at 1:29
• You should have burn-in step, meaning that discard say first 10% of all iterations. You should also check if it is well mixed, for example use ListLinePlot to see this. Accepting ratio should be around %50. You can tune this by choosing different length of interval for UniDist. Is there any reason to choose UniDist? NormalDist is the way to go usually. Feb 21, 2018 at 5:42

OK so here is my understanding of the situation. Be aware that it's the first time I heard about Metropolis or Metropolis-Hasting algorithms: I'm not an expert.

First, a point an terminology: the distinction between Metropolis and Metropolis-Hasting does not seem very clear; for example there is no Wikipedia article on Metropolis algorithm alone, and in this pdf what you call Metropolis is called Metropolis-Hasting.

Let's first start with the implementation of the detailed algorithm you provided for Metropolis:

f[x_] = 0.8 Exp[-0.8 x]/(1 - Exp[-0.8 Pi]);
w[x_] = UnitBox[x/Pi - .5]*f[x];
g[x0_] := RandomVariate[NormalDistribution[x0, 1]]
x0 = RandomReal[{0, Pi}];
stepM[x0_] := Block[{},
xtrial = g[x0];
r = w[xtrial]/w[x0];
eta = RandomReal[{0, 1}];
If[eta <= r, x1 = xtrial;, x1 = x0 ];
x1]

Show[Histogram[NestList[stepM, x0, 100000], 50, "PDF"], Plot[w[x], {x, 0, Pi}]]


Now, let's implement Metropolis-Hasting:

g[x0_] := RandomVariate[NormalDistribution[x0, 1]]
alpha[x0_, xtrial_] := Min[1, w[xtrial]/w[x0]*g[x0]/g[xtrial]]
stepMH[x0_] := Block[{},
xtrial = g[x0];
u = RandomReal[{0, 1}];
If[u <= alpha[x0, xtrial], x1 = xtrial, x1 = x0];
x1]

Show[Histogram[NestList[stepMH, x0, 100000], 50, "PDF"],


Plot[w[x], {x, 0, Pi}]]

As you can see, it's really not as good. And if I take a uniform distribution for g (g[x0_] := RandomVariate[UniformDistribution[{x0 - h/2, x0 + h/2}]]), with h=5 (purposely chosen large), the result is even worse:

What that shows, is that the result is sensitive to the proposal distribution g. The above-linked document recommend using a normal distribution, with a variance chosen has the hessian of w. I'm not sure what it means in this context, but with $\sigma=1.5$ the result is not too bad.

The question that remains open for me is why does MH not use instead alpha[x0_, xtrial_] := Min[1, w[xtrial]/w[x0]] (without the *g[x0]/g[xtrial]), that works better and gives better results. That would be worth asking on stats.stackexchange.com.

• Your code is simply metropolis not metropolis hastings. Look at the step 4 of my two code blocks and you will see the difference. Feb 22, 2018 at 21:36
• @MathDude Could you add the algo for step 4 so that I can make an independent implementation? The question is about Metropolis-Hasting but you provide another algo, which already works for you. Feb 22, 2018 at 21:42
• Here's the new step 4. r = $\dfrac{w(xtrial)}{w(x0)} \dfrac{g(x0|xtrial)}{g(xtrial|x0)}$ where g is the proposal distribution. Feb 23, 2018 at 3:39