plotters_statistical/stats/
normal.rs

1//! Standard-normal inverse CDF (quantile function), for theoretical quantiles
2//! on Q–Q plots.
3
4/// Inverse standard-normal CDF (probit) via Acklam's rational approximation,
5/// accurate to roughly `1e-9` across `(0, 1)`.
6///
7/// Returns `-inf` at `p <= 0` and `+inf` at `p >= 1`.
8///
9/// ```
10/// # use plotters_statistical::stats::norm_ppf;
11/// assert!((norm_ppf(0.5)).abs() < 1e-9);
12/// assert!((norm_ppf(0.975) - 1.959963985).abs() < 1e-6);
13/// ```
14// The coefficients are Acklam's published constants; keep them verbatim.
15#[allow(clippy::excessive_precision)]
16pub fn norm_ppf(p: f64) -> f64 {
17    const A: [f64; 6] = [
18        -3.969683028665376e+01,
19        2.209460984245205e+02,
20        -2.759285104469687e+02,
21        1.383577518672690e+02,
22        -3.066479806614716e+01,
23        2.506628277459239e+00,
24    ];
25    const B: [f64; 5] = [
26        -5.447609879822406e+01,
27        1.615858368580409e+02,
28        -1.556989798598866e+02,
29        6.680131188771972e+01,
30        -1.328068155288572e+01,
31    ];
32    const C: [f64; 6] = [
33        -7.784894002430293e-03,
34        -3.223964580411365e-01,
35        -2.400758277161838e+00,
36        -2.549732539343734e+00,
37        4.374664141464968e+00,
38        2.938163982698783e+00,
39    ];
40    const D: [f64; 4] = [
41        7.784695709041462e-03,
42        3.224671290700398e-01,
43        2.445134137142996e+00,
44        3.754408661907416e+00,
45    ];
46    const P_LOW: f64 = 0.02425;
47    let p_high = 1.0 - P_LOW;
48
49    if p <= 0.0 {
50        return f64::NEG_INFINITY;
51    }
52    if p >= 1.0 {
53        return f64::INFINITY;
54    }
55    if p < P_LOW {
56        let q = (-2.0 * p.ln()).sqrt();
57        (((((C[0] * q + C[1]) * q + C[2]) * q + C[3]) * q + C[4]) * q + C[5])
58            / ((((D[0] * q + D[1]) * q + D[2]) * q + D[3]) * q + 1.0)
59    } else if p <= p_high {
60        let q = p - 0.5;
61        let r = q * q;
62        (((((A[0] * r + A[1]) * r + A[2]) * r + A[3]) * r + A[4]) * r + A[5]) * q
63            / (((((B[0] * r + B[1]) * r + B[2]) * r + B[3]) * r + B[4]) * r + 1.0)
64    } else {
65        let q = (-2.0 * (1.0 - p).ln()).sqrt();
66        -(((((C[0] * q + C[1]) * q + C[2]) * q + C[3]) * q + C[4]) * q + C[5])
67            / ((((D[0] * q + D[1]) * q + D[2]) * q + D[3]) * q + 1.0)
68    }
69}
plotters_statistical/stats/normal.rs

plotters_statistical/stats/
normal.rs
