Hi,
The current icdf in torch.distributions.Normal seems to be numerically unstable compared to scipy.stats.norm.
I implemented a torch version by translating the scipy source code to python and it works the same as the scipy version. I wonder if we could integrate the new code to the current torch.distributions.Normal. Thank you!
Could you explain the numerical instability a bit? I.e. where is it observed and how scipy deals with it?
I’m sure a PR would be welcome, if the instability is reproducible or breaks certain use cases.
Numerical instability issues occur when the input to the inverse CDF is very close to 0 or 1.
scipy handles these cases with approximations and here is the screenshot of the comments in scipy’s source code:
One example is when y=1-1e-8:
torch.distributions.Normal(0, 1).icdf(torch.tensor(y)) returns tensor(inf) (because torch.erfinv returns tensor(inf));scipy.stats.norm.ppf(y) returns 5.612001243305505.Thanks for the update!
Would you like to create a feature request on GitHub and post the explanation there, so that we could discuss it there with the module owners and track it?
Btw. thanks for digging into this. 
that’s kinda related to float32 precision, as 1e-8 < 2**-23 = 1.19e-7
and
torch.erfinv(torch.tensor([1-2**-23])) → 3.7439
issue is this:
np.float32(1)-np.float32(1e-8)
1.0
np.frexp(np.float32(1)-np.float32(1e-8))
(0.5, 1)
I see, this is great, thank you @googlebot !!
Hi @ptrblck, I think there’s no need for a feature request since using DoubleTensor fixes the issue.
Thank you both for your help and suggestions~!
DoubleTensor does not solve the issue with erfinv; it fails at 1e-18. But torch.special.ndtri works to 1e-40 with float32, and to <1e-100 with float64.