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.