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$.
rinvGaussfromSuppDistsorrinvgaussfromstatmoddirectly? $\endgroup$