Accuracy of normal_lccdf much worse than normal_lcdf

[lciti]

Description

The implementation for normal_lcdf has been improved significantly following #1284 by introducing a polynomial approximation of the Mills ratio based on a paper by Cody (1969). It looks like normal_lccdf has not been updated accordingly and instead still uses the older and much less accurate solution.

Incidentally, I also noticed that normal_lcdf and std_normal_lcdf implement the same solution with very similar code. I wonder if one could retain the current implementation of std_normal_lcdf and write std_normal_lccdf, normal_lcdf and normal_lccdf to simply call std_normal_lcdf with appropriate parameters (-x, (x-mu)/sigma, (mu-x)/sigma, respectively).

Current Version:

v3.2.0

[bbbales2]

@PhilClemson what do you think?

[PhilClemson]

Sorry, this is something that was mentioned in the original pull request for the normal_lcdf changes but seems it was forgotten about.

@nhuurre mentioned that the functions could be standardised in a log_Phi function (#1411 (comment)) but this would require separate numerical approximations.

It seems std_normal_lcdf was added by @mcol to follow the structure of the other probability distributions. I guess this might be a wider issue since there's also the std_normal sampling statement within stan?

In any case I don't think there's anything to stop us from adding the improvements to normal_lccdf so that should probably be done as a minimum.

[bob-carpenter]

call std_normal_lcdf with appropriate parameters

What parameters? As a distribution, std-normal doesn't have any parameters. It's just

std_normal_lpdf(y) == normal_lpdf(y | 0, 1).

Accuracy of normal_lccdf much worse than normal_lcdf

lccdf should be more accurate when the cdf is near 1, and the lcdf should be more accurate when the cdf is near 0.

Most of our lcdf and lccdf functions are terrible log(cdf) and log(1 - cdf) implementations, where the lccdf form will be terrible when the cdf is near 1 rather than better.

[lciti]

@bob-carpenter Sorry, I should have been clearer. I meant "parameters" in a CS sense as in arguments. There could be a common C++ function (which @nhuurre called log_Phi) implementing the logCDF of a std normal and all the 4 relevant Stan functions could call that same C++ function under the hood, with appropriate function arguments. For example std_normal_lccdf could call this function with argument -x.

This is a log cdf so I am not sure your argument about the accuracy near 0 or 1 really applies. The domain is not (0,1) but (-inf,0). We do not need to sacrifice the accuracy on either side. The issue with CDF and CCDF is that when they are near 1 we have that the important digits are pushed away on the right of the binary representation because of the "big" (implicit) one on the left. This is not the case for logCDF or logCCDF as when the CDF is close to zero then logCDF and logCCDF become large and negative (retaining all their significant digits if computed right as in @PhilClemson 's implementation); on the other side, when the CDF is near one, logCDF and logCCDF are near zero and can retain all their significant digits if implemented right (as they do in @PhilClemson 's implementation). If we did log(CDF(x)) or log(1-CDF(x)), then it would be a problem.

The current implementation of std_normal_lcdf(x) works extremely well. I am currently writing my Stan code using std_normal_lcdf(-x) (which is no big deal) but I would normally write it as std_normal_lccdf(x), which is simpler (I also have to leave a comment for myself to explain why I am doing the former rather than the latter so I don't risk changing it in the future).

[lciti]

@PhilClemson

@nhuurre mentioned that the functions could be standardised in a log_Phi function (#1411 (comment)) but this would require separate numerical approximations.

I think the comment "this would require separate numerical approximations" only applies when you try to extend log_Phi to non-normal distributions (e.g. exp_mod_normal_lpdf etc). Within the four functions std_normal_lcdf, std_normal_lccdf, normal_lcdf and normal_lccdf, the same approximation works fine as these are just affine transformations of the argument.

[bob-carpenter]

I meant "parameters" in a CS sense as in arguments.

This gets confusing in stats. But whatever you want to call it, std_normal_lpdf is a function from real numbers to real numbers.

The move from std_normal to normal is the usual location/scale generalization. That's also used for most of our other univariate distributions, though we don't have std_logistic or std_student_t.

The whole point of having std_normal was to have a more efficient way of doing things that didn't involve subtracting 0 and dividing by 1. So presumably we're talking about converting normal to std_normal by (y - mu) / sigma, not the other way around. That will get the right density for the normal, but the derivatives are the bigger concern. We don't want to just autodiff (y - mu) / sigma.

I am not sure your argument about the accuracy near 0 or 1 really applies.

As the cdf approaches 0, log(cdf) approaches -infinity and log(1 - cdf) approaches 0. So you're right that we should be able to get decent accuracy both ways. I was tied up in our current implementations, which are a mess.

Edit: I fat-fingered submitting this before finishing.

We do need the complementary forms for the arguments, though. Just consider a beta(y | 0.5, 0.5) distribution. That's going to go to 0 or 1 as y goes to 0 or 1, but you can't represent a y very close to 1 in floating point, so you have to use the lcdf or cdf as appropriate just for the arguments. Same with any distribution truncated at 0, which needs the cdf.

So I don't think we can just drop the lcdf and lccdf distinction computationally.

[bob-carpenter]

The normal is symmetric around the location parameter, so we can relate

normal_lcdf(y | mu, sigma) == normal_lccdf(mu + -(y - mu) | mu, sigma).

[lciti]

@bob-carpenter

So presumably we're talking about converting normal to std_normal by (y - mu) / sigma, not the other way around.

Yes.

That will get the right density for the normal, but the derivatives are the bigger concern. We don't want to just autodiff (y - mu) / sigma.

I see.

We do need the complementary forms for the arguments, though. Just consider a beta(y | 0.5, 0.5) distribution. That's going to go to 0 or 1 as y goes to 0 or 1, but you can't represent a y very close to 1 in floating point, so you have to use the lcdf or cdf as appropriate just for the arguments. Same with any distribution truncated at 0, which needs the cdf.

I was only referring to the Gaussian case. This is particularly useful for me because I can model an exponential distribution as a Gaussian in unconstrained space and then define my constrained variable just as y=std_normal_lcdf(-x). It works much better in my case (in terms of convergence) than using the default log transformation you get with <lower=0>.

[jaburgoyne]

This issue has been dormant for some time, but given the importance of the normal distribution, it might be worth re-prioritising. Especially during initialisation or in some survival modelling, the sampler can end up far out into the tails of a normal distribution. There is a trivial workaround by using log_1m(normal_lcdf(...)), but one needs to have been around the Stan community long enough to know to use it.

Given the time that has passed, it is probably worth having a look across normal_lcdf(), std_normal_lcdf(), normal_lccdf(), and `std_normal_lccdf() to make sure that they are all consistent with each other before closing this issue (see also #2470 and potentially #1025).

Depending on where that lands, the documentation will then need to be updated to reflect what the final version of what these functions actually do and their support regions (see also stan-dev/docs#352).

I'm happy to help if I can: it's just not entirely clear from the history what the actual consensus is/was surrounding which function should be the heart of this family and which implementation should be preferred.

[bob-carpenter]

If there had been a consensus on what to do, it would be reflected here. We'd be happy to take a PR to make any of these functions more numerically accurate, however it's done, as long as the resulting code is human-understandable.