An underflow problem in slice sampling

In a MCMC seminar that Peter and I run, we were recently discussing slice sampling. Here is an overview from Radford Neal himself (open access). Let me skip the details of the algorithm, and pose the particular problem we encountered.

Problem: Compute the energy \(U^*\) of a system. We don't have any idea of the magnitude of this energy, or its sign. Suppose the probability of the system to be in that state is \(\exp(-U^*)\).

This number is going to be extremely small if \(U^*\) is large (which is probable, if you don't work hard at estimating a good initial guess).

The slice sampling algorithm asks us to draw a number \[y \sim U(0, \exp(-U^*)),\] and compute the difference \[\Delta U_y = -\log(y) - U^*.\]
Straightforward Implementation:

The following Octave program computes 10,000 samples of \(\Delta U_y\) to produce the histogram

%
% Small Ustar = 10, no problems
%
Us  = 10;
num = 10000;
y   = rand(num,1) * exp(-Us);
Uy  = -log(y);
dUy = Uy - Us;
[hx, x] = hist(dUy,25);
hx = hx/num;
bar(x, hx)

Histrogram of dUy with U*=10 is an exponential distribution

Unfortunately if I set Us = 1000 in the program above, it fails. The failure is simple to parse. If Us is large, exp(-Us) is below "realmin" in double precision, and is approximated to zero, which spawns a cascade of disastrous events, ending in the logarithm of zero.

There are multiple ways of fixing this problem.

Solution 1:

Since we are really interested in \(\Delta U_y\), we can draw it directly from the correct distribution. In this case, a little analysis suggests that the distribution of \(\Delta U_y\) is an exponential distribution with parameter = 1. Hence, we can skip the auxillary variables \(y\) and \(U_y\) and go directly to the meat of the problem:

Us  = 1000;

num = 10000;

dUy = exprnd(1,num,1);

Solution 2:

Since we observe that \(U^*\) doesn't participate in the distribution, we can think of another method that does not directly involve sampling from an exponential random number generator. We define \(y' \sim U(0,1),\) and \(y = y' * \exp(-U^*)\).

We can write \(-\log(y) = -\log(y') + U^*\).

This implies  that  \(\Delta U_y = -\log(y) - U^* = -\log y'\).

Thus,

yp  = rand(num,1);

dUy = -log(yp);

also works well. The advantage is that one doesn't need an special exponential random number generator.

On my computer, solution 2 was more efficient than solution 1 for num < 50000, after which solution 1 became more efficient.