aprender-core 0.60.0

Next-generation machine learning library in pure Rust
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107


/// Continued fraction for incomplete beta (Lentz's algorithm).
fn beta_continued_fraction(a: f32, b: f32, x: f32) -> f32 {
    let max_iter = 100;
    let eps = 1e-7;

    let qab = a + b;
    let qap = a + 1.0;
    let qam = a - 1.0;

    let mut c = 1.0;
    let mut d = 1.0 - qab * x / qap;
    if d.abs() < 1e-30 {
        d = 1e-30;
    }
    d = 1.0 / d;
    let mut h = d;

    for m in 1..=max_iter {
        let m_f = m as f32;
        let m2 = 2.0 * m_f;

        // Even step
        let aa = m_f * (b - m_f) * x / ((qam + m2) * (a + m2));
        d = 1.0 + aa * d;
        if d.abs() < 1e-30 {
            d = 1e-30;
        }
        c = 1.0 + aa / c;
        if c.abs() < 1e-30 {
            c = 1e-30;
        }
        d = 1.0 / d;
        h *= d * c;

        // Odd step
        let aa = -(a + m_f) * (qab + m_f) * x / ((a + m2) * (qap + m2));
        d = 1.0 + aa * d;
        if d.abs() < 1e-30 {
            d = 1e-30;
        }
        c = 1.0 + aa / c;
        if c.abs() < 1e-30 {
            c = 1e-30;
        }
        d = 1.0 / d;
        let del = d * c;
        h *= del;

        if (del - 1.0).abs() < eps {
            break;
        }
    }

    h
}

/// Log-gamma function `ln Γ(z)` (Lanczos approximation).
///
/// PMAT-904: the raw-space `gamma()` below evaluates `(...).exp()` directly,
/// which overflows f32 (Γ(36) ≈ 1e40, Inf at z ≳ 36). p-value paths for
/// degrees-of-freedom df ≳ 72 must combine gamma terms in LOG space (divide
/// gammas as a subtraction of `ln_gamma`, then a single final `exp`) so no
/// intermediate ever overflows. `ln Γ` grows only ~ z ln z, never overflowing
/// f32 for any df reachable here.
fn ln_gamma(z: f32) -> f32 {
    if z < 0.5 {
        // Reflection: ln Γ(z) = ln(π / sin(πz)) − ln Γ(1−z).
        let sin_pz = (PI * z).sin();
        (PI / sin_pz).ln() - ln_gamma(1.0 - z)
    } else {
        // Lanczos series, evaluated entirely in log space.
        let z = z - 1.0;
        let tmp = z + 5.5;
        let tmp = (z + 0.5) * tmp.ln() - tmp;
        let ser = 1.0 + 76.180_09_f32 / (z + 1.0) - 86.505_32_f32 / (z + 2.0)
            + 24.014_1_f32 / (z + 3.0)
            - 1.231_739_5_f32 / (z + 4.0)
            + 0.001_208_58_f32 / (z + 5.0)
            - 0.000_005_363_82_f32 / (z + 6.0);
        tmp + (ser * (2.0 * PI).sqrt()).ln()
    }
}

/// Gamma function approximation (Stirling's approximation).
///
/// NOTE: this raw-space form overflows f32 for z ≳ 36. Production p-value code
/// uses `ln_gamma` and combines in log space (PMAT-904); this raw form is kept
/// only for the small-argument unit tests.
#[cfg(test)]
fn gamma(z: f32) -> f32 {
    if z < 0.5 {
        // Reflection formula: Γ(z) = π / (sin(πz) * Γ(1-z))
        PI / ((PI * z).sin() * gamma(1.0 - z))
    } else {
        // Stirling's approximation
        let z = z - 1.0;
        let tmp = z + 5.5;
        let tmp = (z + 0.5) * tmp.ln() - tmp;
        let ser = 1.0 + 76.180_09_f32 / (z + 1.0) - 86.505_32_f32 / (z + 2.0)
            + 24.014_1_f32 / (z + 3.0)
            - 1.231_739_5_f32 / (z + 4.0)
            + 0.001_208_58_f32 / (z + 5.0)
            - 0.000_005_363_82_f32 / (z + 6.0);
        (tmp + ser.ln()).exp() * (2.0 * PI).sqrt()
    }
}