# Independence-Metropolis-Hastings Algorithm

IMHA is an importance-sampling version of MCMC, where the proposal is drawn from a fixed distribution g. Usually, g is chosen to be an approximation to f. The acceptance probability becomes

$r(x,y)=\min\{\frac{f(y)g(x)}{f(x)g(y)},1\}$.

A random variable $Z$ has an inverse Gaussian distribution if it has density:

$f(z) \propto z^{-3/2}exp \{ -\theta_1z-\frac{\theta_2}{z} + 2 \sqrt{\theta_1 \theta_2} + log \sqrt{2\theta_2}\}, z>0$,

where $\theta_1 >0$ and $\theta_2 >0$ are parameters. It can be shown that

$E(Z) = \sqrt {\frac{\theta_2}{\theta_1}}$ and $E(\frac{1}{Z})= \sqrt{{\frac{\theta_1}{\theta_2}} + \frac{1}{2 \theta_2}}$.

Say I let $\theta_1 =1.5$ and $\theta_2 =2$. I want to draw a sample of size 1,000 using the independent IMHA. I want to use Gamma distribution as proposal density.

How could I accomplish this using R? And would it be appropriate to assess the accuracy by comparing the means of $Z$ and $1/Z$ from the sample to the theoretical means.

• Any reason not to use rinvGauss from SuppDists or rinvgauss from statmod directly? Mar 26 '15 at 8:51
• Thank you, @ChristophHanck. I had never seen that package before. Are you suggesting something like: rinvgauss(1000, c(1.5,2))? Would that be a correct way to approach this, and if not, how could I use rinvgauss. Mar 26 '15 at 10:04

This is Example 7.18 (pp.293-294) in Monte Carlo Statistical Methods (2004, second edition). Is this a coincidence?!

In any case, here is the whole explanation verbatim from the book:

The inverse Gaussian distribution has the density $$f(z|\theta_1,\theta_2) \propto z^{-3/2} \; \exp \left\{- \theta_1z - {\theta_2 \over z} + 2 \sqrt{\theta_1\theta_2} + \log \; \sqrt{2\theta_2} \right\}$$ on ${\mathbb{R}}_{+}$ $(\theta_1>0,\theta_2>0)$. Denoting $\psi(\theta_1,\theta_2) = 2 \sqrt{\theta_1\theta_2} + \log \sqrt{2\theta_2}$, it follows from a classical result on exponential families (see Brown, 1986, Chapter 2, or Robert, 2001, Lemma 3.3.7) that \begin{eqnarray*} \displaystyle{\mathbb{E}\left[\left(Z,1 / Z\right)\right]} &=& \nabla \psi(\theta_1,\theta_2) \\ &=& \displaystyle{\left(\sqrt{{\theta_2 \over \theta_1}}, \sqrt{{\theta_1 \over \theta_2}} + {1 \over 2\theta_2}\right)} \;. \end{eqnarray*} A possible choice for the simulation of this distribution is the Gamma distribution ${\mathcal G}a(\alpha,\beta)$ in the Metropolis--Hastings algorithm, taking $\alpha = \beta \sqrt{\theta_2/\theta_1}$ so that the means of both distributions coincide. Since $${f(x) \over g(x)} \propto x^{-\alpha-1/2} \; \exp \left\{(\beta-\theta_1)x - {\theta_2 \over x}\right\} \;,$$ the ratio $f/g$ is bounded for $\beta < \theta_1$. The value of $x$ which maximizes the ratio is the solution of $$(\beta-\theta_1)x^2 - \left(\alpha+{1 \over 2}\right)x + \theta_2 = 0 \;;$$ that is, $$x^{*}_\beta = {(\alpha+1/2) - \sqrt{(\alpha+1/2)^2 + 4\theta_2(\theta_1-\beta)} \over 2(\beta-\theta_1)} \;.$$ The analytical optimization (in $\beta$) of $$M(\beta) = (x^{*}_\beta)^{-\alpha-1/2} \; \exp \left\{(\beta-\theta_1)x^{*}_\beta - {\theta_2 \over x^{*}_\beta}\right\}$$ is not possible, although, in this specific case the curve $M(\beta)$ can be plotted for given values of $\theta_1$ and $\theta_2$ and the optimal value $\beta^\star$ can be approximated numerically. Typically, the influence of the choice of $\beta$ must be assessed empirically; that is, by approximating the acceptance rate $\rho$ via the method described above.

$$\begin{matrix} \hline \beta & 0.2 & 0.5 & 0.8 & 0.9 & 1 & 1.1 & 1.2 & 1.5 \cr \hline {\hat\rho(\beta)} &0.22 &0.41 &0.54 &0.56 &0.60 &0.63 &0.64 &0.71 \cr \mathbb{E}[Z] &1.137 &1.158 &1.164 &1.154 &1.133 &1.148 &1.181 &1.148 \cr \mathbb{E}[1/Z] &1.116 &1.108 &1.116 &1.115 &1.120 &1.126 &1.095 &1.115 \cr \hline \end{matrix}$$

Estimation of the means of $Z$ and of $1/Z$ for the inverse Gaussian distribution $\mathcal{IN} (\theta_1,\theta_2)$ by the Metropolis--Hastings algorithm and evaluation of the acceptance rate for the instrumental distribution ${\mathcal G}a(\sqrt{\theta_2 / \theta_1} \; > \beta,\beta)$ \ $(\theta_1=1.5, \theta_2=2$, and $m=5000)$.

Note that a new sample $(y_1,\ldots,y_m)$ must be simulated for every new value of $\beta$. Whereas $y \sim {\mathcal G}a(\alpha,\beta)$ is equivalent to $\beta y \sim {\mathcal G}a(\alpha,1)$, the factor $\alpha$ depends on $\beta$ and it is not possible to use the same sample for several values of $\beta$. The table above provides an evaluation of the rate $\rho$ as a function of $\beta$ and gives estimates of the means of $Z$ and $1/Z$ for $\theta_1 = 1.5$ and $\theta_2 = 2$. The constraint on the ratio $f/g$ then imposes $\beta<1.5$. The corresponding theoretical values are respectively $1.155$ and $1.116$, and the optimal value of $\beta$ is $\beta^*=1.5$.

• May I take this post to imply that you have access to the tex file of the book :-)? Mar 26 '15 at 10:07
• @ChristophHanck: you may. Mar 26 '15 at 10:09