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

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?).

• All you need to do is divide each gamma variate by their sum. How, then, does underflow occur? And how does taking the logarithm solve this problem (you can't compute the sum without exponentiating back again anyway)?
– whuber
Commented Mar 7, 2011 at 19:54
• @whuber In log space, you compute the sum and then subtract it from each element. So, this avoids the first point of underflow. There is a bit of further processing when these dirichlets serve as mixture components and get multiplied by small numbers again. Commented Mar 7, 2011 at 19:58
• Adding the logs is mathematically incorrect: it corresponds to multiplying the gammas rather than adding them. Yes, you might get working results, but they definitely will not have a Dirichlet distribution! Again, what exactly is the nature of the original underflow and what calculations are you doing when it happens? What are the actual values you are working with?
– whuber
Commented Mar 7, 2011 at 20:02
• @whuber I might have simplified a bit too much in my description. I do forall i { t = gamma(a,b); sum += t; d[i] = log(t)}; logsum = log(sum); forall i { d[i] -= logsum; }. Previously, this underflowed if a was very small. Commented Mar 7, 2011 at 20:05
• Got it: for $\alpha$ near 0 you're going to be in trouble no matter what. Interesting problem!
– whuber
Commented Mar 7, 2011 at 21:15

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.

• Just an argument to make your edit a bit more elegant, you don't really need to appeal to integration here. Just use the fact that $\frac{\Gamma(\alpha)}{\Gamma(\alpha)+\Gamma(1)}\sim Beta(\alpha,1)$, plus that $\Gamma(\alpha)+\Gamma(1)\sim \Gamma(\alpha+1)$. These are both standard properties of the beta and gamma distributions. Also, when $\alpha\approx 0$ we have roughly $y\sim expo(1)$, which may be quicker to simulate ($-\log(u)$) than a general $\Gamma(\alpha+1)$ random variable. Commented Oct 9, 2011 at 0:42

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.

• Could you explain what gsl_rng_uniform_pos and gsl_ran_gamma do? I would guess the first returns a uniform random value between 0 and r and the second is related to a Gamma(1+a, b) value--maybe it's an incomplete Gamma? Overall this looks very close to the approximation I suggested (except, in reviewing it, it's obvious I forgot to specify the divide by $\alpha$ part, which is essential!)
– whuber
Commented Mar 15, 2011 at 22:59
• I edited my answer to include more detail now. Commented Mar 16, 2011 at 18:15
• Thank you: but what is "r"? (Note that the recursion is bounded: at most one recursive call will be made, because a > 0 implies 1.0+a > 1.)
– whuber
Commented Mar 16, 2011 at 18:35
• r is the random number generator (where you are getting the random numbers from). Commented Mar 16, 2011 at 19:11
• Ah, this is clever: the product of a $\Gamma(\alpha+1)$ and an independent $B(\alpha,1)$ variate turns out to be a $\Gamma(\alpha)$ variate. I edited my reply so it points to your solution and explains why it works.
– whuber
Commented Mar 16, 2011 at 20:58