statrs-dev / statrs

Statistical computation library for Rust
https://docs.rs/statrs/latest/statrs/
MIT License
582 stars 83 forks source link

normal cdf function gives probability of 0.0 for very low probabilities #191

Closed robinpaul85 closed 1 month ago

robinpaul85 commented 1 year ago

I have noticed that for very low probability the normal cdf function gives a probability of 0.0 when actual probability are in the order of 10^(-30)

WarrenWeckesser commented 1 year ago

I can't reproduce this. If I evaluate the CDF of a normal distribution with mean 0 and standard deviation 1 at x = -35.0, I get 1.124910706417e-268.

Can you provide code that gives a minimal working example that shows the unexpected result?

robinpaul85 commented 1 year ago

let z = 8.78491962525629; let n = Normal::new(0.0, 1.0).unwrap(); let p_g = 1.0 - n.cdf(z);

This gives p_g = 0.0

In comparison, when I use the r_stats crate, documentation here I get:

let p_g2 = r_stats::normal_cdf(z, 0.0, 1.0, false, false);

I get p_g2 = 0.0000000000000000007823818993898402

WarrenWeckesser commented 1 year ago

The expression 1 - cdf(x) is often called the complementary CDF or the survival function. For values of x that are sufficiently large, actually computing it with the expression 1 - cdf(x) is not a good method, as you have observed. The problem is that with x that large, cdf(x) is numerically 1.0. So you end up with 1 - cdf(x) = 1 - 1 = 0.

Fortunately the sf method was added to the continuous distributions in version 0.16 of statrs, so if you are using the latest version, you can compute the desired quantity with the new method. Using the names from your code sample, that would be p_g = n.sf(z);.

If you can't upgrade statrs, a simple fix is to use the symmetry of the normal distribution: sf(μ + x) = cdf(μ - x).