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
98
99
100
101
102
103

// Translation of nmath's dnorm
//
// 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, Positive64, Real64};

use crate::traits::DPQ;

/// Compute the density of the normal distribution.
pub fn dnorm<R1: Into<Real64>, R2: Into<Real64>, P: Into<Positive64>>(
    x: R1,
    mu: R2,
    sigma: P,
    log: bool,
) -> Real64 {
    let mut x = x.into().unwrap();
    let mu = mu.into().unwrap();
    let sigma = sigma.into().unwrap();

    if !sigma.is_finite() {
        return f64::d_0(log).into();
    }
    if !x.is_finite() && mu == x {
        return f64::nan().into();
    }
    if sigma == 0.0 {
        /* sigma == 0 */
        return if x == mu {
            f64::infinity()
        } else {
            f64::d_0(log)
        }
        .into();
    }
    x = (x - mu) / sigma;

    if !x.is_finite() {
        return f64::d_0(log).into();
    }

    x = x.abs();
    if x >= 2.0 * f64::MAX.sqrt() {
        return f64::d_0(log).into();
    }
    if log {
        return (-(strafe_consts::LN_SQRT_2TPI + 0.5 * x * x + sigma.ln())).into();
    }
    // M_1_SQRT_2PI = 1 / sqrt(2 * pi)
    // more accurate, less fast :
    if x < 5.0 {
        return (strafe_consts::_1DSQRT_2TPI * (-0.5 * x * x).exp() / sigma).into();
    } else {
        /* x*x  may lose upto about two digits accuracy for "large" x
        * Morten Welinder's proposal for PR#15620
        * https://bugs.r-project.org/show_bug.cgi?id=15620

        * -- 1 --  No hoop jumping when we underflow to zero anyway:

        *  -x^2/2 <         log(2)*.Machine$double.min.exp  <==>
        *     x   > sqrt(-2*log(2)*.Machine$double.min.exp) =IEEE= 37.64031
        * but "thanks" to denormalized numbers, underflow happens a bit later,
        *  effective.D.MIN.EXP <- with(.Machine, double.min.exp + double.ulp.digits)
        * for IEEE, DBL_MIN_EXP is -1022 but "effective" is -1074
        * ==> boundary = sqrt(-2*log(2)*(.Machine$double.min.exp + .Machine$double.ulp.digits))
        *              =IEEE=  38.58601
        * [on one x86_64 platform, effective boundary a bit lower: 38.56804]
        */
        if x > (-2.0
            * std::f64::consts::LN_2
            * (f64::MIN_EXP as f64 + 1.0 - std::f64::MANTISSA_DIGITS as f64))
            .sqrt()
        {
            return 0.0.into();
        }

        /* Now, to get full accuracy, split x into two parts,
        *  x = x1+x2, such that |x2| <= 2^-16.
        * Assuming that we are using IEEE doubles, that means that
        * x1*x1 is error free for x<f64::MAX_EXP (but we have x < 38.6 anyway).

        * If we do not have IEEE this is still an improvement over the naive formula.
        */
        let x1 = x.ldexp(16.0).round().ldexp(-16.0);
        let x2 = x - x1;
        let ret = strafe_consts::_1DSQRT_2TPI / sigma
            * ((-0.5 * x1 * x1).exp() * ((-0.5 * x2 - x1) * x2).exp());
        ret.into()
    }
}