0
$\begingroup$

The target distribution is of the form:

$ p(x) = x^{-6}.e^{\frac{-2.475}{x}}$ with a support in the interval $[0.0, 2.0]$.

This gives a plot like

plot target over [0.0, 2.0]

Now, to choose a proposal kernel, I think a lognormal may suffice. To determine what the values of the hyperparameters for the lognormal is that could be a best fit to the target distribution, I played around and identify s = 0.5, loc = 0, scale = 0.5 and a constant = 0.30 that multiples the proposal kernel so that it looks like the target distribution.

So my proposal kernel has the form $q(.|.)$ = 0.3 * lognorm.pdf(x, s, loc, scale) using the package https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.lognorm.html

Comparing the target distribution with the introduced proposal kernel $q(.|.) $ (this proposal kernel will now be used for the Metropolis - Hasting sampling):

enter image description here

Here is how I implement the M-H algorithm:

  1. First I initialise a value for x (call this $x_{0}$)
  2. I specify a markov chain length $t \in [1, T]$
  3. Then I condition the proposal kernel $q(.|.)$ on $x_{0}$ to sample an updated value for $x$: this is $x^{t} \sim q(x^{t} | x_{0})$. This is reflected at the algorithmic level by replacing the loc parameter in $q(.|.)$ with $x_{0}$
  4. Next, sample $x_{t-1}$ from $q(x^{t-1} | x_{t})$. Algorithmically, this is replacing the loc parameter in $q(.|.)$ with $x_{t}$
  5. 3., 4. allows me to compute the Hasting ratio $H$ for $x$: this is just $x_{t-1}/x{t}$
  6. the target distribution $p$ of $x_{t}$ and $p$ of $x_{t-1}$ respectively: $p(x_{t}), p(x_{t-1})$
  7. the target distribution ratio $T$: $p(x_{t})/p(x_{t-1})$
  8. compute the M-H ratio $r$: $r = T * H$
  9. randomly sample $u$ from a uniform distribution on the support $(0,1)$
  10. if $u >= r$ then accept $x_{t}$ else reject.
  11. The next iteration for $t$ begins: this time $x_{t}$ is used as an initial value in 3.

After 1000 t's the plot below emerge:

enter image description here

matplotlib.pyplot.hist(accepted_x_samples, density = False, weights = weights)

As you can see the shape looks OK but the bin axis is off. The vertical axis is also suspect.

Been working on this for days diligently but I cannot figure this out. Please, any help is appreciated.

$\endgroup$

1 Answer 1

1
$\begingroup$

There are several issues here:

  • There are various versions of the Metropolis--Hastings algorithm. Some draw samples "around" the current point (in your notation, this corresponds to changing loc to $x_{t-1}$), others have a fixed proposal distribution that does not change with your current location $x_{t}$ (fixed loc). Both work! However, since you fitted your proposal to be close to your target, it seems that the second approach is what you want.
  • Please distinguish between a sample $x$ from a some distribution with density $\rho$ and the density value $\rho(x)$ at that point. In 5. you divide $x_{t-1}/x_{t}$, this is not meaningful. You have to divide the corresponding proposal density values!
  • In 4. you sample $x_{t-1}$. This is not how Metropolis--Hastings works: You use the point $x_{t-1}$ that you already have from the previous step ($x_{0}$ in the first step) and evaluate the proposal densities correspondingly.
  • (You didn't write out what "reject" means: You set $x_{t}$ to equal $x_{t-1}$.)
  • Since you fitted your distribution so nicely, importance sampling is something you might consider, depending on what you want to do with your samples.

Hope this helps!

$\endgroup$
2
  • $\begingroup$ Thank you for the suggestions. To your second point, that was a typo - I divided the conditional density in my code. With MH sampling, what motivates a good proposal density if not a density that resembles the target? $\endgroup$ Feb 22 at 15:42
  • $\begingroup$ If you have a density that resembles the target, then it is typically a good choice, but use it without changing loc to $x_{t-1}$. In your particular example, the proposal is supported on $\mathbb{R}^{+}$. If you change loc to $x_{t-1}$ in each step, you can only move to the right and never to the left! This explains your observations. I suggest that you have a look at your output sequence. Given your construction, it should be strictly increasing, which is definitely not what you want for a meaningful sampling method, right? $\endgroup$ Feb 22 at 23:48

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .