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

// Translation of nmath's dt
//
// 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 strafe_type::{FloatConstraint, Rational64, Real64};

use crate::{
    distribution::{
        func::{bd0, stirlerr},
        norm::dnorm,
    },
    traits::DPQ,
};

/// The t density is evaluated as
/// (n/2).sqrt() / ((n+1)/2) * Gamma((n+3)/2) / Gamma((n+2)/2).
/// * (1+x^2/n)^(-n/2)
/// / \sqrt( 2 pi (1+x^2/n) )
///
/// This form leads to a stable computation for all
/// values of n, including n -> 0 and n -> infinity.
pub fn dt<RE: Into<Real64>, RA: Into<Rational64>>(x: RE, n: RA, log: bool) -> Real64 {
    let x = x.into().unwrap();
    let n = n.into().unwrap();

    // := ((1 + x2n).sqrt()).ln() = (1 + x2n).ln()/2
    if !x.is_finite() {
        return f64::d_0(log).into();
    }
    if !n.is_finite() {
        return dnorm(x, 0.0, 1.0, log);
    }

    let mut u = 0.0;
    let t = -bd0(n / 2.0, (n + 1.0) / 2.0).unwrap() + stirlerr((n + 1.0) / 2.0) - stirlerr(n / 2.0);
    let x2n = x * x / n;
    let mut ax = 0.0;
    let mut l_x2n = 0.0;
    let lrg_x2n = x2n > 1.0 / 2.2204460492503131e-16;
    if lrg_x2n {
        // large x^2/n :
        ax = x.abs(); // = x2n.ln()/2 = 1/2 * (x^2 / n).ln()
        l_x2n = ax.ln() - n.ln() / 2.0;
        u = n * l_x2n
    } else if x2n > 0.2 {
        l_x2n = (1.0 + x2n).ln() / 2.0;
        u = n * l_x2n
    } else {
        l_x2n = x2n.ln_1p() / 2.0;
        u = -bd0(n / 2.0, (n + x * x) / 2.0).unwrap() + x * x / 2.0
    }

    //old: return  R_D_fexp(M_2PI*(1+x2n), t-u);
    // R_D_fexp(f,x) :=  (give_log ? -0.5*f.ln()+(x) : x.exp()/f.sqrt())
    // f = 2pi*(1+x2n)
    // ==> 0.5*f.ln() = 2pi.ln()/2 + (1+x2n).ln()/2 = 2pi.ln()/2 + l_x2n
    // 1/f.sqrt() = 1/(2pi * (1+ x^2 / n)).sqrt()
    // = 1/2pi.sqrt()/(|x|/n.sqrt()*(1+1/x2n).sqrt())
    // = M_1_SQRT_2PI * n.sqrt()/ (|x|*(1+1/x2n).sqrt())
    if log {
        return (t - u - (strafe_consts::LN_SQRT_2TPI + l_x2n)).into();
    }

    // else :  if(lrg_x2n) : (1 + 1/x2n).sqrt() ='= 1.sqrt() = 1
    let I_sqrt_ = if lrg_x2n {
        (n.sqrt()) / ax
    } else {
        (-l_x2n).exp()
    };
    return ((t - u).exp() * strafe_consts::_1DSQRT_2TPI * I_sqrt_).into();
}