How to quickly sample X if exp(X) ~ Gamma?

Question

I have a simple sampling problem, where my inner loop looks like:

v = sample_gamma(k, a)

where sample_gamma samples from the Gamma distribution to form a Dirichlet sample.

It works well, but for some values of k/a, some of the downstream computation underflows.

I adapted it to use log space variables:

v = log(sample_gamma(k, a))

After adapting all of the rest of the programme, it works correctly (at least it gives me the same exact results on test cases). However, it is slower than before.

Is there a way to directly sample $X, \exp(X) \sim \text{Gamma}$ without using slow functions like $\log()$? I tried googling for this, but I do not even know if this distribution has a common name (log-gamma?).

Answer 1

Consider a small shape parameter $\alpha$ near 0, such as $\alpha = 1/100$. In the range between 0 and $\alpha$, $e^{-\alpha}$ is approximately $1$, so the Gamma pdf is approximately $x^{\alpha-1}dx / \Gamma(\alpha)$. This can be integrated to an approximate CDF, $F_\alpha(x) = \frac{x^\alpha}{\alpha \Gamma(\alpha)}$. Inverting it, we see a $1/\alpha$ power: a huge exponent. For $\alpha = 1/100$ this causes some chance of underflow (a double precision value less than $10^{-300}$, more or less). Here is a plot of the chance of getting underflow as a function of the base-ten logarithm of $\alpha$:

One solution is to exploit this approximation for generating log(Gamma) variates: in effect, try to generate a Gamma variate and if it's too small, generate its logarithm from this approximate power distribution (as shown below). (Do this repeatedly until the log is within the underflow range, so that it is a valid substitute for the original underflowing variate.) For the Dirichlet calculation, subtract the maximum of all the logarithms from each of the log values: this implicitly rescales all the Gamma variates so it won't affect the Dirichlet values. Treat any resulting log that is too small (say, less than -100) as being the log of a true zero. Exponentiate the other logs. Now you can proceed without underflow.

This is going to take even longer than before, but at least it will work!

To generate an approximate log Gamma variate with shape parameter $\alpha$, precompute $C = \log(\Gamma(\alpha)) + \log(\alpha)$. This is easy, because there are algorithms to compute values of log Gamma directly. Generate a uniform random float between 0 and 1, take its logarithm, divide by $\alpha$, and add $C$ to it.

Because the scale parameter merely rescales the variate, there is no problem accommodating it in these procedures. You don't even need it if all scale parameters are the same.

Edit

In another reply the OP describes a method in which the $1/\alpha$ power of a uniform variate (a $B(\alpha)$ variate) is multiplied by a $\Gamma(\alpha+1)$ variate. This works because the pdf of the joint distribution of these two variates equals $\left(\alpha x^{\alpha-1}\right) \left(y^{\alpha}e^{-y}dy/\Gamma(\alpha+1)\right)$. To find the pdf of $z = xy$ we substitute $y \to z/x$, divide by the Jacobean $x$, and integrate out $x$. The integral must range from $z$ to $\infty$ because $0 \le y \le 1$, whence

$$\text{pdf}(z)=\frac{\alpha}{\Gamma(\alpha+1)}\int_z^{\infty}{\left(x^\alpha / x\right) e^{-x} (z/x)^{\alpha-1} dx} dz = \frac{1}{\Gamma(\alpha)}z^{\alpha-1}e^{-z}dz,$$

which is the pdf of a $\Gamma(\alpha)$ distribution.

The whole point is that when $0 \lt \alpha \lt 1$, a value drawn from $\Gamma(\alpha+1)$ is unlikely to underflow and by summing its log and $1/\alpha$ times the log of an independent uniform variate we will have the log of a $\Gamma(\alpha)$ variate. The log is likely to be very negative, but we will have bypassed the construction of its antilog, which will underflow in a floating point representation.

Answer 2

I'm answering my own question, but I found a pretty good solution, even if I do not fully understand it. Looking at the code from the GNU Scientific Library, here is how it samples gamma variables (r is the random number generator, a is $\alpha$ and b is $\beta$):

  if (a < 1)
    {
      double u = gsl_rng_uniform_pos (r);
      return gsl_ran_gamma (r, 1.0 + a, b) * pow (u, 1.0 / a);
   }

gsl_ran_gamma is the function which returns a gamma random sample (so the above is a recursive call), while gsl_rng_uniform_pos returns a uniformly distributed number in $(0,1)$ (the _pos is for strictly positive as it is guaranteed to not return 0.0).

Therefore, I can take the log of the last expression and use

return log(gsl_ran_gamma(r, 1.0 + a, b)) + log(u)/a;

To get what I wanted. I now have two log() calls (but one less pow()), but the result is probably better. Before, as whuber pointed out, I had something raised to the power of $1/a$, potentially a huge number. Now, in logspace, I'm multiplying by $1/a$. So, it is less likely to underflow.