Repository navigation
Accuracy of normal_lccdf much worse than normal_lcdf #1985
Description
Activity
@PhilClemson what do you think?
Sorry, this is something that was mentioned in the original pull request for the
normal_lcdfchanges but seems it was forgotten about.@nhuurre mentioned that the functions could be standardised in a
log_Phifunction (#1411 (comment)) but this would require separate numerical approximations.It seems
std_normal_lcdfwas added by @mcol to follow the structure of the other probability distributions. I guess this might be a wider issue since there's also thestd_normalsampling statement within stan?In any case I don't think there's anything to stop us from adding the improvements to
normal_lccdfso that should probably be done as a minimum.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.
@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 examplestd_normal_lccdfcould 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))orlog(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 usingstd_normal_lcdf(-x)(which is no big deal) but I would normally write it asstd_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).@nhuurre mentioned that the functions could be standardised in a
log_Phifunction (#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_Phito non-normal distributions (e.g.exp_mod_normal_lpdfetc). Within the four functionsstd_normal_lcdf,std_normal_lccdf,normal_lcdfandnormal_lccdf, the same approximation works fine as these are just affine transformations of the argument.I meant "parameters" in a CS sense as in arguments.
This gets confusing in stats. But whatever you want to call it,
std_normal_lpdfis a function from real numbers to real numbers.The move from
std_normaltonormalis the usual location/scale generalization. That's also used for most of our other univariate distributions, though we don't havestd_logisticorstd_student_t.The whole point of having
std_normalwas 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 convertingnormaltostd_normalby(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.
The normal is symmetric around the location parameter, so we can relate
normal_lcdf(y | mu, sigma) == normal_lccdf(mu + -(y - mu) | mu, sigma).Reacted by Luca CitiSo 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 defaultlogtransformation you get with<lower=0>.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.
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.
- linked a pull request that will close this issueAligned normal_lccdf and std_normal_lccdf with the lcdf functions #3363
on Aug 29, 2026 Thanks so much for adding, @jaburgoyne, and welcome to the Stan project!
Further thanks to everyone who ran the hackathon at StanCon that helped out.
Description
The implementation for
normal_lcdfhas been improved significantly following #1284 by introducing a polynomial approximation of the Mills ratio based on a paper by Cody (1969). It looks likenormal_lccdfhas not been updated accordingly and instead still uses the older and much less accurate solution.Incidentally, I also noticed that
normal_lcdfandstd_normal_lcdfimplement the same solution with very similar code. I wonder if one could retain the current implementation ofstd_normal_lcdfand writestd_normal_lccdf,normal_lcdfandnormal_lccdfto simply callstd_normal_lcdfwith appropriate parameters (-x, (x-mu)/sigma, (mu-x)/sigma, respectively).Current Version:
v3.2.0