cas_compute/funcs/
probability.rs

1//! Probability density and distribution functions.
2
3use cas_attrs::builtin;
4use crate::consts::{E, ONE, PI, TAU, TWO, ZERO};
5use crate::primitive::{float, float_from_str, int};
6use once_cell::sync::Lazy;
7use rug::{float::Special, ops::Pow, Float, Integer};
8use super::miscellaneous::partial_factorial;
9
10/// The error function, `erf(x)`.
11#[derive(Debug)]
12pub struct Erf;
13
14#[cfg_attr(feature = "numerical", builtin)]
15impl Erf {
16    pub fn eval_static(x: Float) -> Float {
17        x.erf()
18    }
19}
20
21/// The complementary error function, `erfc(x)`.
22#[derive(Debug)]
23pub struct Erfc;
24
25#[cfg_attr(feature = "numerical", builtin)]
26impl Erfc {
27    pub fn eval_static(x: Float) -> Float {
28        x.erfc()
29    }
30}
31
32/// Precomputed values for the inverse error function.
33static INV_ERF_A: Lazy<[Float; 4]> = Lazy::new(|| [
34    float_from_str("0.886226899"),
35    float_from_str("-1.645349621"),
36    float_from_str("0.914624893"),
37    float_from_str("-0.140543331"),
38]);
39
40static INV_ERF_B: Lazy<[Float; 5]> = Lazy::new(|| [
41    float(&*ONE),
42    float_from_str("-2.118377725"),
43    float_from_str("1.442710462"),
44    float_from_str("-0.329097515"),
45    float_from_str("0.012229801"),
46]);
47
48static INV_ERF_C: Lazy<[Float; 4]> = Lazy::new(|| [
49    float_from_str("-1.970840454"),
50    float_from_str("-1.62490649"),
51    float_from_str("3.429567803"),
52    float_from_str("1.641345311"),
53]);
54
55static INV_ERF_D: Lazy<[Float; 3]> = Lazy::new(|| [
56    float(&*ONE),
57    float_from_str("3.543889200"),
58    float_from_str("1.637067800"),
59]);
60
61/// The inverse error function, `inverf(x)`.
62#[derive(Debug)]
63pub struct Inverf;
64
65#[cfg_attr(feature = "numerical", builtin)]
66impl Inverf {
67    pub fn eval_static(x: Float) -> Float {
68        if x <= -1 {
69            return float(Special::NegInfinity);
70        } else if x >= 1 {
71            return float(Special::Infinity);
72        }
73
74        // implementation from http://www.mimirgames.com/articles/programming/approximations-of-the-inverse-error-function/
75        let mut z = float(x.abs_ref());
76
77        let fold_arr = |arr: &[Float], factor: &Float| {
78            (1..arr.len())
79                .rev()
80                .fold(float(&arr[arr.len() - 1]), |acc, i| {
81                    acc * factor + float(&arr[i - 1])
82                })
83        };
84
85        let mut r = if z <= 0.7 {
86            let x2 = float(z.square_ref());
87            let num = &z * fold_arr(&*INV_ERF_A, &x2);
88            let den = fold_arr(&*INV_ERF_B, &x2);
89            num / den
90        } else {
91            let y = {
92                let inner_log = (&*ONE - float(&z)) / &*TWO;
93                (-inner_log.ln()).sqrt()
94            };
95            let num = fold_arr(&*INV_ERF_C, &y);
96            let den = fold_arr(&*INV_ERF_D, &y);
97            num / den
98        };
99
100        r = if x.is_sign_negative() { -r } else { r };
101        z = if x.is_sign_negative() { -z } else { z };
102
103        // double newton's method
104        let f = 2 / float(&*PI).sqrt();
105        r -= (float(r.erf_ref()) - &z) / (&f * float((-float(r.square_ref())).exp_ref()));
106        r -= (float(r.erf_ref()) - &z) / (&f * float((-float(r.square_ref())).exp_ref()));
107
108        r
109    }
110}
111
112/// Normal probability density function.
113///
114/// This function does not actually represent any probability. To get the probability of a normally
115/// distributed random variable being between a value `a` and `b`, use [`Normcdf`].
116#[derive(Debug)]
117pub struct Normpdf;
118
119#[cfg_attr(feature = "numerical", builtin)]
120impl Normpdf {
121    pub fn eval_static(x: Float, m: Option<Float>, s: Option<Float>) -> Float {
122        let m = m.unwrap_or_else(|| float(&*ZERO));
123        let s = s.unwrap_or_else(|| float(&*ONE));
124        let exp_arg = Float::exp(float(x - m).square() / (-2 * float(s.square_ref())));
125        let bot = s * float(TAU.sqrt_ref());
126        exp_arg / bot
127    }
128}
129
130/// Normal cumulative distribution function.
131///
132/// Returns the probability that a random variable, chosen from a normal distribution with mean `m`
133/// and standard deviation `s`, will be between `a` and `b`, inclusive.
134#[derive(Debug)]
135pub struct Normcdf;
136
137#[cfg_attr(feature = "numerical", builtin)]
138impl Normcdf {
139    pub fn eval_static(a: Float, b: Float, m: Option<Float>, s: Option<Float>) -> Float {
140        let m = m.unwrap_or_else(|| float(&*ZERO));
141        let s = s.unwrap_or_else(|| float(&*ONE));
142        let sqrt_two = float(TWO.sqrt_ref());
143        let z_a = (a - &m) / float(&s * &sqrt_two);
144        let z_b = (b - &m) / float(&s * &sqrt_two);
145        (z_b.erf() - z_a.erf()) / &*TWO
146    }
147}
148
149/// Inverse normal cumulative distribution function.
150///
151/// Returns the value `x` such that `normcdf(-∞, x, m, s) = p`, i.e., the value `x` such that the
152/// probability of a random variable, chosen from a normal distribution with mean `m` and standard
153/// deviation `s`, being between `-∞` and `x` is `p`.
154#[derive(Debug)]
155pub struct Invnorm;
156
157#[cfg_attr(feature = "numerical", builtin)]
158impl Invnorm {
159    pub fn eval_static(p: Float, m: Option<Float>, s: Option<Float>) -> Float {
160        let m = m.unwrap_or_else(|| float(&*ZERO));
161        let s = s.unwrap_or_else(|| float(&*ONE));
162        let sqrt_two = float(TWO.sqrt_ref());
163        let z = sqrt_two * Inverf::eval_static(2 * p - 1);
164        m + s * z
165    }
166}
167
168/// Geometric probability function.
169///
170/// Returns the probability of the first success of an event occurring on the `n`th trial, where the
171/// probability of success on a single trial is `p`.
172#[derive(Debug)]
173pub struct Geompdf;
174
175#[cfg_attr(feature = "numerical", builtin)]
176impl Geompdf {
177    pub fn eval_static(p: Float, n: Integer) -> Float {
178        if n <= *ZERO {
179            return float(&*ZERO);
180        }
181
182        let q = float(&*ONE - &p);
183        p * q.pow(n - 1)
184    }
185}
186
187/// Cummulative geometric probability function.
188///
189/// Returns the probability of the first success of an event occurring on or before the `n`th
190/// trial, where the probability of success on a single trial is `p`.
191#[derive(Debug)]
192pub struct Geomcdf;
193
194#[cfg_attr(feature = "numerical", builtin)]
195impl Geomcdf {
196    pub fn eval_static(p: Float, n: Integer) -> Float {
197        if n <= *ZERO {
198            return float(&*ZERO);
199        }
200
201        let q = float(&*ONE - &p);
202        &*ONE - q.pow(n)
203    }
204}
205
206/// Binomial probability function.
207///
208/// Returns the probability of exactly `x` successes occurring in `n` trials, where the probability
209/// of success on a single trial is `p`.
210#[derive(Debug)]
211pub struct Binompdf;
212
213#[cfg_attr(feature = "numerical", builtin)]
214impl Binompdf {
215    pub fn eval_static(n: Integer, p: Float, x: Integer) -> Float {
216        if x < *ZERO || x > n {
217            return float(&*ZERO);
218        }
219
220        let q = float(&*ONE - &p);
221
222        let c = q.pow(n.clone() - &x);
223        let b = p.pow(&x);
224        let a = super::combinatoric::Ncr::eval_static(n, x);
225        a * b * c
226    }
227}
228
229/// Cummulative binomial probability function.
230///
231/// Returns the probability of `x` or fewer successes occurring in `n` trials, where the
232/// probability of success on a single trial is `p`.
233#[derive(Debug)]
234pub struct Binomcdf;
235
236#[cfg_attr(feature = "numerical", builtin)]
237impl Binomcdf {
238    pub fn eval_static(n: Integer, p: Float, mut x: Integer) -> Float {
239        if x < *ZERO {
240            return float(&*ZERO);
241        } else if x >= n {
242            return float(&*ONE);
243        }
244
245        let mut sum = float(&*ZERO);
246        while x >= *ZERO {
247            sum += Binompdf::eval_static(n.clone(), p.clone(), x.clone());
248            x -= 1;
249        }
250        sum
251    }
252}
253
254/// Poisson probability function.
255///
256/// Returns the probability of exactly `k` events occurring in an arbitrary time interval, where
257/// the mean number of occurrences of the event in the time interval is `l`.
258#[derive(Debug)]
259pub struct Poisspdf;
260
261#[cfg_attr(feature = "numerical", builtin)]
262impl Poisspdf {
263    pub fn eval_static(k: Integer, l: Float) -> Float {
264        if k < *ZERO {
265            return float(&*ZERO);
266        }
267
268        let b = float(&*E).pow(&*l.as_neg());
269        let a = l.pow(&k);
270        let c = partial_factorial(k, int(1));
271        a * b / c
272    }
273}
274
275/// Cummulative poisson probability function.
276///
277/// Returns the probability of `k` or fewer events occurring in an arbitrary time interval, where
278/// the mean number of occurrences of the event in the time interval is `l`.
279#[derive(Debug)]
280pub struct Poisscdf;
281
282#[cfg_attr(feature = "numerical", builtin)]
283impl Poisscdf {
284    pub fn eval_static(k: Integer, l: Float) -> Float {
285        if k < *ZERO {
286            return float(&*ZERO);
287        }
288
289        fn reg_gamma(s: Float, x: Float) -> Float {
290            float(s.gamma_inc_ref(&x)) / s.gamma()
291        }
292
293        reg_gamma(float(k + 1), l)
294    }
295}
cas_compute/funcs/probability.rs

cas_compute/funcs/
probability.rs
