r2rs-nmath 0.1.1

Statistics programming for Rust based on R's nmath package
Documentation 
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

// Translation of nmath's df
//
// Copyright (C) 2020 neek-sss
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program.  If not, see <https://www.gnu.org/licenses/>.

use num_traits::Float;
use strafe_type::{FloatConstraint, Rational64, Real64};

use crate::{
    distribution::{binom::dbinom_raw, gamma::dgamma},
    traits::DPQ,
};

/// The density function of the F distribution.
/// To evaluate it, write it as a Binomial probability with p = x*m/(n+x*m).
/// For m >= 2, we use the simplest conversion.
/// For m < 2, (m-2)/2 < 0 so the conversion will not work, and we must use
/// a second conversion.
/// Note the division by p; this seems unavoidable
/// for m < 2, since the F density has a singularity as x (or p) -> 0.
pub fn df<RE: Into<Real64>, RA1: Into<Rational64>, RA2: Into<Rational64>>(
    x: RE,
    m: RA1,
    n: RA2,
    log: bool,
) -> Real64 {
    let x = x.into().unwrap();
    let m = m.into().unwrap();
    let n = n.into().unwrap();

    let mut p = 0.0;
    let mut q = 0.0;
    let mut f = 0.0;
    let mut dens = 0.0;

    if x < 0.0 {
        return f64::d_0(log).into();
    }
    if x == 0.0 {
        return if m > 2.0 {
            f64::d_0(log)
        } else if m == 2.0 {
            f64::d_1(log)
        } else {
            f64::infinity()
        }
        .into();
    }
    if !m.is_finite() && !n.is_finite() {
        /* both +Inf */
        return if x == 1.0 {
            f64::infinity()
        } else {
            f64::d_0(log)
        }
        .into();
    }
    if !n.is_finite() {
        /* must be +Inf by now */
        return dgamma(x, m / 2.0, 2.0 / m, log);
    }
    if m > 1e14 {
        /* includes +Inf: code below is inaccurate there */
        dens = dgamma(1.0 / x, n / 2.0, 2.0 / n, log).unwrap();
        return if log {
            dens - 2.0 * x.ln()
        } else {
            dens / (x * x)
        }
        .into();
    }

    f = 1.0 / (n + x * m);
    q = n * f;
    p = x * m * f;

    if m >= 2.0 {
        f = m * q / 2.0;
        dens = dbinom_raw((m - 2.0) / 2.0, (m + n - 2.0) / 2.0, p, q, log)
    } else {
        f = m * m * q / (2.0 * p * (m + n));
        dens = dbinom_raw(m / 2.0, (m + n) / 2.0, p, q, log)
    }
    if log { (f.ln()) + dens } else { f * dens }.into()
}