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
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):
Here is how I implement the M-H algorithm:
- First I initialise a value for x (call this $x_{0}$)
- I specify a markov chain length $t \in [1, T]$
- 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}$
- Next, sample $x_{t-1}$ from $q(x^{t-1} | x_{t})$. Algorithmically, this is replacing the loc parameter in $q(.|.)$ with $x_{t}$
- 3., 4. allows me to compute the Hasting ratio $H$ for $x$: this is just $x_{t-1}/x{t}$
- the target distribution $p$ of $x_{t}$ and $p$ of $x_{t-1}$ respectively: $p(x_{t}), p(x_{t-1})$
- the target distribution ratio $T$: $p(x_{t})/p(x_{t-1})$
- compute the M-H ratio $r$: $r = T * H$
- randomly sample $u$ from a uniform distribution on the support $(0,1)$
- if $u >= r$ then accept $x_{t}$ else reject.
- 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:
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.