statest/
numerical.rs

1//! Numerical calculation for tests.
2
3use num_traits::Float;
4use statrs::distribution::Univariate;
5
6/// Calculate derivative.
7pub fn deriv<S, T>(dist: &T, x: S, h: S) -> S
8where
9    S: Float,
10    T: Univariate<S, S>
11{
12    let fx = dist.cdf(x);
13    deriv_with_fx(dist, x, fx, h)
14}
15
16/// Calculate derivative with fx value which was already calculated.
17pub fn deriv_with_fx<S, T>(dist: &T, x: S, fx: S, h: S) -> S
18where
19    S: Float,
20    T: Univariate<S, S>
21{
22    let fx_h = dist.cdf(x + h);
23    (fx_h - fx) / h
24}
25
26/// Newton method to solve equations.
27pub fn newton<S, T>(dist: &T, y: S, mu: S, epsilon: S) -> S 
28where
29    S: Float,
30    T: Univariate<S, S>
31{
32    let fx0 = dist.cdf(mu);
33    if fx0 + epsilon >= y && fx0 - epsilon <= y {
34        return mu;
35    }
36    let mut x = if fx0 > mu {
37        mu + S::epsilon()
38    } else {
39        mu - S::epsilon()
40    };
41    loop {
42        let fx = dist.cdf(x);
43        if fx + epsilon >= y && fx - epsilon <= y {
44            break;
45        }
46        let df = if x >= mu {
47            deriv_with_fx(dist, x, fx, epsilon)
48        } else {
49            deriv_with_fx(dist, x, fx, -epsilon)
50        };
51        x = x + (y - fx) / df;
52    }
53    x
54}