rsomics_quantile_transform/
ndtri.rs

1//! Inverse normal CDF via a Cephes `ndtri` port.
2//!
3//! `scipy.special.ndtri` is the Cephes `ndtri` (Moshier). QuantileTransformer
4//! uses `stats.norm.ppf` (= `ndtri`) for the normal output distribution, then
5//! clips to `[ndtri(BOUNDS_THRESHOLD - eps), ndtri(1 - (BOUNDS_THRESHOLD - eps))]`
6//! where `BOUNDS_THRESHOLD = 1e-7` and `eps = np.spacing(1)`.
7//!
8//! The central `|y-0.5| ≤ 3/8` branch and the two `z = sqrt(-2 ln y)` tail
9//! branches match scipy bit-for-bit on the same architecture. Cross-arch (arm
10//! vs x86) the last bits of the transcendental can differ by ≤1 ULP; compat
11//! tests use ≤1e-12 relative tolerance for normal output.
12
13// Coefficients verbatim from Cephes at full source precision.
14#![allow(clippy::excessive_precision)]
15
16const S2PI: f64 = 2.506_628_274_631_000_502_42;
17
18const P0: [f64; 5] = [
19    -5.996_335_010_141_078_952_67e1,
20    9.800_107_541_859_996_615_36e1,
21    -5.667_628_574_690_702_934_39e1,
22    1.393_126_093_872_796_795_03e1,
23    -1.239_165_838_673_812_580_16e0,
24];
25const Q0: [f64; 8] = [
26    1.954_488_583_381_417_598_34e0,
27    4.676_279_128_988_815_384_53e0,
28    8.636_024_213_908_905_905_75e1,
29    -2.254_626_878_541_193_705_27e2,
30    2.002_602_123_800_606_603_59e2,
31    -8.203_722_561_683_333_399_12e1,
32    1.590_562_251_262_116_955_15e1,
33    -1.183_316_211_213_300_031_42e0,
34];
35
36const P1: [f64; 9] = [
37    4.055_448_923_059_624_199_23e0,
38    3.152_510_945_998_938_661_54e1,
39    5.716_281_922_464_212_881_62e1,
40    4.408_050_738_932_008_347_00e1,
41    1.468_495_619_288_580_240_14e1,
42    2.186_633_068_507_902_675_39e0,
43    -1.402_560_791_713_544_958_75e-1,
44    -3.504_246_268_278_482_034_18e-2,
45    -8.574_567_851_546_854_136_11e-4,
46];
47const Q1: [f64; 8] = [
48    1.577_998_832_564_667_497_31e1,
49    4.539_076_351_288_792_105_84e1,
50    4.131_720_382_546_720_304_40e1,
51    1.504_253_856_929_075_034_08e1,
52    2.504_649_462_083_094_159_79e0,
53    -1.421_829_228_547_877_885_74e-1,
54    -3.808_064_076_915_782_771_94e-2,
55    -9.332_594_808_954_574_273_72e-4,
56];
57
58const P2: [f64; 9] = [
59    3.237_748_917_769_460_359_70e0,
60    6.915_228_890_689_842_116_95e0,
61    3.938_810_252_924_744_434_15e0,
62    1.333_034_608_158_075_423_89e0,
63    2.014_853_895_491_790_815_38e-1,
64    1.237_166_348_178_200_213_58e-2,
65    3.015_815_535_082_354_160_07e-4,
66    2.658_069_746_867_375_508_32e-6,
67    6.239_745_391_849_832_937_30e-9,
68];
69const Q2: [f64; 8] = [
70    6.024_270_393_647_420_142_55e0,
71    3.679_835_638_561_608_594_03e0,
72    1.377_020_994_890_813_302_71e0,
73    2.162_369_935_944_966_358_90e-1,
74    1.342_040_060_885_431_890_37e-2,
75    3.280_144_646_821_277_391_04e-4,
76    2.892_478_647_453_806_839_36e-6,
77    6.790_194_080_099_812_744_25e-9,
78];
79
80/// Inverse standard normal CDF: the `x` for which Φ(x) = `y`. Cephes `ndtri`.
81#[must_use]
82#[allow(clippy::manual_range_contains)]
83pub fn ndtri(y0: f64) -> f64 {
84    if y0 == 0.0 {
85        return f64::NEG_INFINITY;
86    }
87    if y0 == 1.0 {
88        return f64::INFINITY;
89    }
90    if y0 < 0.0 || y0 > 1.0 {
91        return f64::NAN;
92    }
93
94    let mut code = true;
95    let mut y = y0;
96    if y > 1.0 - 0.135_335_283_236_612_691_89 {
97        y = 1.0 - y;
98        code = false;
99    }
100
101    if y > 0.135_335_283_236_612_691_89 {
102        y -= 0.5;
103        let y2 = y * y;
104        let x = y + y * (y2 * polevl(y2, &P0) / p1evl(y2, &Q0));
105        return x * S2PI;
106    }
107
108    let x = (-2.0 * y.ln()).sqrt();
109    let x0 = x - x.ln() / x;
110    let z = 1.0 / x;
111    let x1 = if x < 8.0 {
112        z * polevl(z, &P1) / p1evl(z, &Q1)
113    } else {
114        z * polevl(z, &P2) / p1evl(z, &Q2)
115    };
116    let x = x0 - x1;
117    if code { -x } else { x }
118}
119
120fn polevl(x: f64, coef: &[f64]) -> f64 {
121    let mut ans = coef[0];
122    for &c in &coef[1..] {
123        ans = ans * x + c;
124    }
125    ans
126}
127
128fn p1evl(x: f64, coef: &[f64]) -> f64 {
129    let mut ans = x + coef[0];
130    for &c in &coef[1..] {
131        ans = ans * x + c;
132    }
133    ans
134}
135
136/// `stats.norm.ppf(BOUNDS_THRESHOLD - np.spacing(1))` — the clip floor for
137/// normal output. Pre-computed constant; cross-arch ndtri of this input is
138/// bit-identical because the argument is deep in a tail, fully polynomial.
139pub const CLIP_MIN: f64 = -5.199_337_582_605_575;
140/// `stats.norm.ppf(1 - (BOUNDS_THRESHOLD - np.spacing(1)))` — the clip ceiling.
141pub const CLIP_MAX: f64 = 5.199_337_582_703_42;
142
143#[cfg(test)]
144mod tests {
145    use super::*;
146
147    fn close(got: f64, want: f64) {
148        let rel = (got - want).abs() / want.abs().max(f64::MIN_POSITIVE);
149        assert!(rel <= 1e-12, "got {got:e} want {want:e} rel={rel:e}");
150    }
151
152    #[test]
153    fn ndtri_known_quantiles() {
154        close(ndtri(0.975), 1.959_963_984_540_054);
155        close(ndtri(0.995), 2.575_829_303_548_900_8);
156        close(ndtri(0.5), 0.0);
157    }
158
159    #[test]
160    fn ndtri_endpoints() {
161        assert_eq!(ndtri(0.0), f64::NEG_INFINITY);
162        assert_eq!(ndtri(1.0), f64::INFINITY);
163        assert!(ndtri(-0.1).is_nan());
164        assert!(ndtri(1.1).is_nan());
165    }
166
167    #[test]
168    fn clip_constants_match_scipy() {
169        // np.spacing(1) == f64::EPSILON (ULP of 1.0 = 2.22e-16)
170        // scipy.stats.norm.ppf(1e-7 - np.spacing(1)) == -5.199337582605575
171        // scipy.stats.norm.ppf(1 - (1e-7 - np.spacing(1))) == 5.19933758270342
172        // Verified via struct.pack('>d', v).hex()
173        let got_min = ndtri(1e-7_f64 - f64::EPSILON);
174        let got_max = ndtri(1.0 - (1e-7_f64 - f64::EPSILON));
175        close(got_min, CLIP_MIN);
176        close(got_max, CLIP_MAX);
177    }
178}
rsomics_quantile_transform/ndtri.rs

rsomics_quantile_transform/
ndtri.rs
