3
$\begingroup$

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.

enter image description here

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

enter image description here

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

enter image description here

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.

enter image description here

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

enter image description here

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}]];
  ]
$\endgroup$
  • $\begingroup$ 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). $\endgroup$ – anderstood Feb 19 '18 at 3:27
  • 1
    $\begingroup$ 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. $\endgroup$ – Szabolcs Feb 19 '18 at 9:41
  • $\begingroup$ You could try reproducing the algorithm in this wikipedia section. $\endgroup$ – anderstood Feb 20 '18 at 20:39
  • $\begingroup$ @Szabolcs I've added comments above which should appear after peer review. Please let me know if it helps. $\endgroup$ – ChemDude Feb 21 '18 at 1:29
  • $\begingroup$ 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. $\endgroup$ – OkkesDulgerci Feb 21 '18 at 5:42
3
+25
$\begingroup$

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

enter image description here

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

enter image description here

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:

enter image description here

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.

$\endgroup$
  • $\begingroup$ Your code is simply metropolis not metropolis hastings. Look at the step 4 of my two code blocks and you will see the difference. $\endgroup$ – MathDude Feb 22 '18 at 21:36
  • $\begingroup$ @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. $\endgroup$ – anderstood Feb 22 '18 at 21:42
  • $\begingroup$ 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. $\endgroup$ – ChemDude Feb 23 '18 at 3:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.