Numerical precision of normal_lcdf

[alecksphillips]

Description

The normal_lcdf(x | mu, sigma) function suffers from precision problems at x >> mu/sigma when compared to the equivalent pnorm(x, mu, sigma, log.p=T) in base R.

Example

R code:

library(rstan)

stan_model_code <- "
functions {
  real test_normal_lcdf(real x) {
    return normal_lcdf(x | 0, 1);
  }
}
parameters{}
model{}"

expose_stan_functions(stanc(model_code = stan_model_code))

x <- seq(7.5,7.51,0.00001)
y_stan <- sapply(x, test_normal_lcdf)
y_R <- pnorm(x,log.p=T)

df <- data.frame(
  x = c(x,x),
  y = c(y_R, y_stan),
  method = factor(c(rep_len("R", length(x)), rep_len("Stan", length(x))))
)

plot(df$x, df$y, col = df$method, xlab = "x", ylab = "normal_lcdf(x)")
legend(x = 7.502, y = -3e-14, legend = levels(df$method), col = c(1:2), pch = 16)

Expected Output:

See plot comparing Stan and R's output here: https://discourse.mc-stan.org/t/numerical-precision-of-normal-lcdf/9685

Current Version:

v2.19.1

[bob-carpenter]

Thanks for reporting.

This is a known problem with pretty much every one of our lcdf and lccdf functions, which are implemented as log(foo_cdf(...)) and log1m(foo_cdf(...)). What's needed are iterative algorithms for the cdfs that work on the log scale with stable derivatives. Often what happens when you get out into the tails, log(foo_cdf(...)) remains finite, but one or more of its derivatives, which have the term 1 / foo_cdf(...), explode.

[avehtari]

I don't think derivatives are the problem here, but computing log(1+erf(z)), using erf from boost.

R stats package pnorm.c https://github.com/SurajGupta/r-source/blob/a28e609e72ed7c47f6ddfbb86c85279a0750f0b7/src/nmath/pnorm.c has numerically more accurate case for computing lcdf in tails with trivial derivatives. both pnorm.c and Boost erf are using the same rational Chebyshev approximation, but pnorm.c handles the cdf and lcdf far in tails with more care (see lines 199-204 and 219-260). Boost has cdf for normal but not lcdf, so it seems that the solution would be the follow the approach used in pnorm.c

[avehtari]

I tested that using log(pnorm(x,log.p=F)) in R results also to loss of precision (with same discrete levels in lcdf values, but slightly different location of discontinuities. So it really is the case that lcdf specific computation is needed (and the pnorm.c has the reference implementation)

[wds15]

How does the phi_approx compare to this? It may not be as precise, but it should be better as it stays continuous...of course i don’t expect great precision in the far tails, but maybe that’s better than the discrete steps you see. To be clear, I mean the common logisitic approximation to the cdf which is described in the stan manuals.

[avehtari]

log(Phi_approx(x)) underflows to 0 in this range, ie, it's much worse.

[wds15]

No...you should use the log inv logit function and define it accordingly.

[avehtari]

You wrote phi_approx and I tested phi_approx. Can you be more specific and how your suggestion compares to the algorithm used in pnorm.c?

[wds15]

You should always avoid calculating on the natural scale and log at the very end if you target log.

So this is phi approx:

inv_logit(0.07056 * pow(x, 3.0) + 1.5976 * x)

And in stan we have a log inv logit...so you should try

Log_inv_logit(0.07056 * pow(x, 3.0) + 1.5976 * x)

I hope this helps.

[avehtari]

I know about "avoid calculating on the natural scale and log at the very end if you target log." and that's why I recommend the approach from pnorm.c which does things in log scale and has been demonstrated to produce good results.

So this is phi approx:

The confusion came as there is also a function called Phi_approx().
I know this approximation, but I doubt this is not good approximation for normal_lcdf compared to what pnorm.c does. I think that in this issue, we should focus on accurate computation for normal_lcdf and continue in the original discourse thread discussion about whether @alecksphillips could use this approximation instead of normal_lcdf.

Phi_approx often works because often there is no difference even between logit or probit (Phi_approx has x^3 term so it's not just logit).

[wds15]

Sure...the issue should be on the pnorm.c implementation and my note is better made in discourse. For practical purposes the log inv logit solution may give good enough results. If the approximation gives already good fits then this issue can be prioritized adequatley.

[avehtari]

Log_inv_logit(0.07056 * pow(x, 3.0) + 1.5976 * x)

Relative error compared to solution in pnorm.c in the example range is about 4550000%

[wds15]

Sure...but this is far in the tails! Does it really matter?

[avehtari]

I think it matters for this issue, and it's separate thing whether it matters in a certain model, which should be discussed somewhere else.