Robust Computation of Phi(b) - Phi(a)

Developers
1

Every once in a while I find I need to compute the log probability of an interval under a normal distribution. Although in principle this is straightforward –

log(P([a,b] | mu, sigma)) = log(Phi((b-mu)/sigma) - Phi((a-mu)/sigma))

– in practice this gets more complex if you want something that is very robust. For example, (b-a)/sigma may be very small, or both a and b may have large magnitudes with the same sign. You’d like to be able to get a high-precision answer whenever that answer is representable as a double. Preferably low relative error, but at least low absolute error (recall that we’re computing the log of the interval probability).

I haven’t seen any special function in Stan for doing this. Is anyone aware of anything in the research literature addressing this problem? I haven’t found anything yet. If I have to I’ll work out my own solution, but I’d rather not re-invent the wheel.

1 Like
2

Is log_diff_exp(normal_lcdf(b | mu, sigma), normal_lcdf(a | mu, sigma) insufficiently robust?

1 Like
3

For a, b > 0 complementary CDF difference is more robust

log_diff_exp(normal_lccdf(a | mu, sigma), normal_lccdf(b | mu, sigma)

When a,b are very close to each other (less than 10^{-5} or so) you might use a Taylor approximation

real x = 0.5*(a+b)/sigma;
real d = a-b;
return log(2/pi()) - log(sigma) - 0.5*square(x) + log(d) + (square(x)-1)*square(d)/6;

but I think log_diff_exp is still more accurate when \left|a\right|,\left|b\right|>\sigma.

Technically, there’s ordered probit distribution but that’s implemented as log_diff_exp.

ordered_probit_lpmf(2| 0, [a,b]) // log(Phi(b) - Phi(a))
6 Likes