ggplot_rs/stat/
distribution.rs

1//! Self-contained quantile / density functions for the theoretical
2//! distributions of `stat_qq` / `stat_qq_band` (normal, Student t,
3//! exponential, half-normal).
4//!
5//! These must work in every feature configuration — including the minimal,
6//! plotters-free `--no-default-features` build used by the DuckDB extensions —
7//! so they don't depend on the optional `regression` feature's distribution
8//! library. Accuracy (checked against R in the tests):
9//!
10//! - `qnorm`: Wichura's AS 241 (PPND16), ~1e-16 relative;
11//! - `pt` / `qt`: regularized incomplete beta (Lentz continued fraction) and a
12//!   safeguarded Newton inversion, ~1e-12;
13//! - `ln_gamma`: Lanczos (g = 7, n = 9), ~1e-15.
14
15// Published algorithm coefficients are kept verbatim.
16#![allow(clippy::excessive_precision)]
17
18use std::f64::consts::PI;
19
20/// Evaluate a polynomial with coefficients in ascending order (Horner).
21fn poly(c: &[f64], x: f64) -> f64 {
22    c.iter().rev().fold(0.0, |acc, &k| acc * x + k)
23}
24
25/// Standard-normal quantile (R's `qnorm(p)`), Wichura's AS 241.
26pub fn qnorm(p: f64) -> f64 {
27    if p.is_nan() || !(0.0..=1.0).contains(&p) {
28        return f64::NAN;
29    }
30    if p == 0.0 {
31        return f64::NEG_INFINITY;
32    }
33    if p == 1.0 {
34        return f64::INFINITY;
35    }
36    const A: [f64; 8] = [
37        3.387_132_872_796_366_608,
38        133.141_667_891_784_377_45,
39        1_971.590_950_306_551_442_7,
40        13_731.693_765_509_461_125,
41        45_921.953_931_549_871_457,
42        67_265.770_927_008_700_853,
43        33_430.575_583_588_128_105,
44        2_509.080_928_730_122_672_7,
45    ];
46    const B: [f64; 8] = [
47        1.0,
48        42.313_330_701_600_911_252,
49        687.187_007_492_057_908_3,
50        5_394.196_021_424_751_107_7,
51        21_213.794_301_586_595_867,
52        39_307.895_800_092_710_61,
53        28_729.085_735_721_942_674,
54        5_226.495_278_852_854_561,
55    ];
56    const C: [f64; 8] = [
57        1.423_437_110_749_683_577_34,
58        4.630_337_846_156_545_295_9,
59        5.769_497_221_460_691_405_5,
60        3.647_848_324_763_204_605_04,
61        1.270_458_252_452_368_382_58,
62        0.241_780_725_177_450_611_77,
63        0.022_723_844_989_269_184_583_3,
64        7.745_450_142_783_414_076_4e-4,
65    ];
66    const D: [f64; 8] = [
67        1.0,
68        2.053_191_626_637_758_821_87,
69        1.676_384_830_183_803_849_4,
70        0.689_767_334_985_100_004_55,
71        0.148_103_976_427_480_074_59,
72        0.015_198_666_563_616_457_196_6,
73        5.475_938_084_995_344_946e-4,
74        1.050_750_071_644_416_843_24e-9,
75    ];
76    const E: [f64; 8] = [
77        6.657_904_643_501_103_777_2,
78        5.463_784_911_164_114_369_9,
79        1.784_826_539_917_291_335_8,
80        0.296_560_571_828_504_891_23,
81        0.026_532_189_526_576_123_093,
82        0.001_242_660_947_388_078_438_6,
83        2.711_555_568_743_487_578_15e-5,
84        2.010_334_399_292_288_132_65e-7,
85    ];
86    const F: [f64; 8] = [
87        1.0,
88        0.599_832_206_555_887_937_69,
89        0.136_929_880_922_735_805_31,
90        0.014_875_361_290_850_614_852_5,
91        7.868_691_311_456_132_591e-4,
92        1.846_318_317_510_054_681_8e-5,
93        1.421_511_758_316_445_888_7e-7,
94        2.044_263_103_389_939_785_64e-15,
95    ];
96    let q = p - 0.5;
97    if q.abs() <= 0.425 {
98        let r = 0.180_625 - q * q;
99        return q * poly(&A, r) / poly(&B, r);
100    }
101    let r = if q < 0.0 { p } else { 1.0 - p };
102    let r = (-r.ln()).sqrt();
103    let v = if r <= 5.0 {
104        let r = r - 1.6;
105        poly(&C, r) / poly(&D, r)
106    } else {
107        let r = r - 5.0;
108        poly(&E, r) / poly(&F, r)
109    };
110    if q < 0.0 {
111        -v
112    } else {
113        v
114    }
115}
116
117/// Standard-normal density.
118pub fn dnorm(x: f64) -> f64 {
119    (-0.5 * x * x).exp() / (2.0 * PI).sqrt()
120}
121
122/// `ln Γ(x)` for `x > 0` (Lanczos approximation, reflection below 0.5).
123pub fn ln_gamma(x: f64) -> f64 {
124    const G: f64 = 7.0;
125    const COEF: [f64; 9] = [
126        0.999_999_999_999_809_93,
127        676.520_368_121_885_1,
128        -1_259.139_216_722_402_8,
129        771.323_428_777_653_13,
130        -176.615_029_162_140_59,
131        12.507_343_278_686_905,
132        -0.138_571_095_265_720_12,
133        9.984_369_578_019_571_6e-6,
134        1.505_632_735_149_311_6e-7,
135    ];
136    if x < 0.5 {
137        // Γ(x)Γ(1−x) = π / sin(πx)
138        return (PI / (PI * x).sin()).abs().ln() - ln_gamma(1.0 - x);
139    }
140    let x = x - 1.0;
141    let mut a = COEF[0];
142    let t = x + G + 0.5;
143    for (i, &c) in COEF.iter().enumerate().skip(1) {
144        a += c / (x + i as f64);
145    }
146    0.5 * (2.0 * PI).ln() + (x + 0.5) * t.ln() - t + a.ln()
147}
148
149/// Continued fraction for the incomplete beta (modified Lentz).
150fn beta_cf(a: f64, b: f64, x: f64) -> f64 {
151    const TINY: f64 = 1e-300;
152    let (qab, qap, qam) = (a + b, a + 1.0, a - 1.0);
153    let mut c = 1.0;
154    let mut d = 1.0 - qab * x / qap;
155    if d.abs() < TINY {
156        d = TINY;
157    }
158    d = 1.0 / d;
159    let mut h = d;
160    for m in 1..=500 {
161        let m = m as f64;
162        let m2 = 2.0 * m;
163        let aa = m * (b - m) * x / ((qam + m2) * (a + m2));
164        d = 1.0 + aa * d;
165        if d.abs() < TINY {
166            d = TINY;
167        }
168        c = 1.0 + aa / c;
169        if c.abs() < TINY {
170            c = TINY;
171        }
172        d = 1.0 / d;
173        h *= d * c;
174        let aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2));
175        d = 1.0 + aa * d;
176        if d.abs() < TINY {
177            d = TINY;
178        }
179        c = 1.0 + aa / c;
180        if c.abs() < TINY {
181            c = TINY;
182        }
183        d = 1.0 / d;
184        let del = d * c;
185        h *= del;
186        if (del - 1.0).abs() < 1e-15 {
187            break;
188        }
189    }
190    h
191}
192
193/// Regularized incomplete beta `I_x(a, b)`.
194pub fn beta_inc(a: f64, b: f64, x: f64) -> f64 {
195    if x <= 0.0 {
196        return 0.0;
197    }
198    if x >= 1.0 {
199        return 1.0;
200    }
201    let ln_front = ln_gamma(a + b) - ln_gamma(a) - ln_gamma(b) + a * x.ln() + b * (1.0 - x).ln();
202    if x < (a + 1.0) / (a + b + 2.0) {
203        ln_front.exp() * beta_cf(a, b, x) / a
204    } else {
205        1.0 - ln_front.exp() * beta_cf(b, a, 1.0 - x) / b
206    }
207}
208
209/// Student-t lower-tail probability `P(T ≤ t)` (R's `pt(t, df)`).
210pub fn pt(t: f64, df: f64) -> f64 {
211    if t.is_nan() || df.is_nan() || df <= 0.0 {
212        return f64::NAN;
213    }
214    if t.is_infinite() {
215        return if t > 0.0 { 1.0 } else { 0.0 };
216    }
217    let x = df / (df + t * t);
218    let tail = 0.5 * beta_inc(df / 2.0, 0.5, x);
219    if t > 0.0 {
220        1.0 - tail
221    } else {
222        tail
223    }
224}
225
226/// Student-t density (R's `dt(x, df)`).
227pub fn dt(x: f64, df: f64) -> f64 {
228    if df <= 0.0 || x.is_nan() {
229        return f64::NAN;
230    }
231    (ln_gamma((df + 1.0) / 2.0)
232        - ln_gamma(df / 2.0)
233        - 0.5 * (df * PI).ln()
234        - (df + 1.0) / 2.0 * (1.0 + x * x / df).ln())
235    .exp()
236}
237
238/// Student-t quantile for any `p` (R's `qt(p, df)`): closed forms for
239/// `df = 1, 2`, otherwise a bracketed Newton inversion of [`pt`].
240pub fn qt(p: f64, df: f64) -> f64 {
241    if p.is_nan() || df.is_nan() || df <= 0.0 || !(0.0..=1.0).contains(&p) {
242        return f64::NAN;
243    }
244    if p == 0.0 {
245        return f64::NEG_INFINITY;
246    }
247    if p == 1.0 {
248        return f64::INFINITY;
249    }
250    if p == 0.5 {
251        return 0.0;
252    }
253    if (df - 1.0).abs() < 1e-12 {
254        return (PI * (p - 0.5)).tan();
255    }
256    if (df - 2.0).abs() < 1e-12 {
257        return (2.0 * p - 1.0) / (2.0 * p * (1.0 - p)).sqrt();
258    }
259    if p > 0.5 {
260        return -qt_lower(1.0 - p, df);
261    }
262    qt_lower(p, df)
263}
264
265/// `qt` for `p < 0.5` (a negative quantile), solved in the lower tail so
266/// tiny `p` keep full precision.
267fn qt_lower(p: f64, df: f64) -> f64 {
268    // Bracket [lo, hi] with pt(lo) < p < pt(hi) = 0.5.
269    let mut hi = 0.0;
270    let mut lo = -1.0;
271    let mut guard = 0;
272    while pt(lo, df) > p && guard < 2000 {
273        hi = lo;
274        lo *= 2.0;
275        guard += 1;
276    }
277    // Start from the normal quantile with a Cornish–Fisher correction.
278    let z = qnorm(p);
279    let g1 = (z.powi(3) + z) / 4.0;
280    let mut t = (z + g1 / df).clamp(lo, hi);
281    for _ in 0..200 {
282        let f = pt(t, df) - p;
283        if f.abs() <= 1e-15 * p.max(1e-300) {
284            break;
285        }
286        if f > 0.0 {
287            hi = t;
288        } else {
289            lo = t;
290        }
291        let d = dt(t, df);
292        let mut next = t - f / d;
293        if !(next.is_finite() && next > lo && next < hi) {
294            next = 0.5 * (lo + hi);
295        }
296        if (next - t).abs() <= 1e-14 * t.abs().max(1.0) {
297            t = next;
298            break;
299        }
300        t = next;
301    }
302    t
303}
304
305#[cfg(test)]
306mod tests {
307    use super::*;
308
309    fn close(a: f64, b: f64, tol: f64) -> bool {
310        (a - b).abs() <= tol * b.abs().max(1.0)
311    }
312
313    #[test]
314    fn qnorm_matches_r() {
315        // R: sprintf("%.17g", qnorm(c(1e-10, 0.001, 0.025, 0.25, 0.6, 0.975, 0.999999)))
316        let cases = [
317            (1e-10, -6.361_340_902_404_056_6),
318            (0.001, -3.090_232_306_167_813_5),
319            (0.025, -1.959_963_984_540_054),
320            (0.25, -0.674_489_750_196_081_74),
321            (0.6, 0.253_347_103_135_799_78),
322            (0.975, 1.959_963_984_540_054),
323            (0.999_999, 4.753_424_308_817_089),
324        ];
325        for (p, want) in cases {
326            assert!(close(qnorm(p), want, 1e-13), "qnorm({p}) = {}", qnorm(p));
327        }
328        assert_eq!(qnorm(0.5), 0.0);
329        assert_eq!(qnorm(0.0), f64::NEG_INFINITY);
330        assert!(qnorm(1.5).is_nan());
331    }
332
333    #[test]
334    fn t_distribution_matches_r() {
335        // R: qt(c(0.01, 0.1, 0.3, 0.975), df = 3); qt(0.001, 7.5); pt(-2.5, 4); dt(1.3, 5)
336        let q3 = [
337            (0.01, -4.540_702_858_568_133),
338            (0.1, -1.637_744_353_696_208_9),
339            (0.3, -0.584_389_727_439_818_65),
340            (0.975, 3.182_446_305_283_707_8),
341        ];
342        for (p, want) in q3 {
343            assert!(
344                close(qt(p, 3.0), want, 1e-10),
345                "qt({p}, 3) = {}",
346                qt(p, 3.0)
347            );
348        }
349        assert!(close(qt(0.001, 7.5), -4.630_319_094_034_099_7, 1e-9));
350        assert!(close(pt(-2.5, 4.0), 0.033_383_272_405_994_063, 1e-11));
351        assert!(close(dt(1.3, 5.0), 0.158_476_735_728_982_41, 1e-12));
352        // Closed forms.
353        assert!(close(qt(0.9, 1.0), 3.077_683_537_175_254_4, 1e-12));
354        assert!(close(qt(0.9, 2.0), 1.885_618_083_164_126_9, 1e-12));
355        // Round trip in the far tail.
356        let t = qt(1e-9, 4.0);
357        assert!(close(pt(t, 4.0), 1e-9, 1e-8), "{t}");
358    }
359
360    #[test]
361    fn ln_gamma_matches_known_values() {
362        assert!(close(ln_gamma(0.5), 0.5 * PI.ln(), 1e-14));
363        assert!(close(ln_gamma(10.0), (362_880.0f64).ln(), 1e-14));
364        assert!(close(ln_gamma(1.0), 0.0, 1e-14));
365    }
366}
ggplot_rs/stat/distribution.rs

ggplot_rs/stat/
distribution.rs
