3
$\begingroup$

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.

$\endgroup$
2
  • $\begingroup$ Any reason not to use rinvGauss from SuppDists or rinvgauss from statmod directly? $\endgroup$ Mar 26 '15 at 8:51
  • $\begingroup$ 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. $\endgroup$
    – user84756
    Mar 26 '15 at 10:04
1
$\begingroup$

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$.


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

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.