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 dpois
//
// 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 2 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, Positive64, Real64};

use crate::{
    distribution::func::{ebd0, stirlerr},
    func::lgamma,
    traits::DPQ,
};

const M_SQRT_2PI: f64 = 2.50662827463100050241576528481104525301; /* sqrt(2*pi) */
const M_2PI: f64 = 6.283185307179586476925286766559;
// sqrt(2 * Rmpfr::Const("pi", 128))
const x_LRG: f64 = 2.86111748575702815380240589208115399625e+307; /* = 2^1023 / pi */

/// dpois_raw() computes the Poisson probability  $ lb^x exp(-lb) / x! $.
/// This does not check that x is an integer, since dgamma() may
/// call this with a fractional x argument. Any necessary argument
/// checks should be done in the calling function.
pub fn dpois_raw(x: f64, lambda: f64, log: bool) -> f64 {
    /*       x >= 0 ; integer for dpois(), but not e.g. for pgamma()!
        lambda >= 0
    */
    if lambda == 0.0 {
        return if x == 0.0 {
            f64::d_1(log)
        } else {
            f64::d_0(log)
        };
    } // including for the case where  x = lambda = +Inf
    if !lambda.is_finite() {
        return f64::d_0(log);
    }
    if x < 0.0 {
        return f64::d_0(log);
    }
    if x <= lambda * 2.2250738585072014e-308 {
        return (-lambda).d_exp(log);
    }
    if lambda < x * 2.2250738585072014e-308 {
        if !x.is_finite() {
            // lambda < x = +Inf
            return f64::d_0(log);
        }
        // else
        return (-lambda + x * lambda.ln() - lgamma(x + 1.0).unwrap()).d_exp(log);
    }
    // R <= 4.0.x  had   return (_2TPI * x).d_fexp(-stirlerr(x) - tof64!(bd0(fin64!(x), fin64!(lambda))), log);
    let mut yh = 0.0;
    let mut yl = 0.0;
    ebd0(x, lambda, &mut yh, &mut yl);
    yl += stirlerr(x);
    let Lrg_x = x >= x_LRG; //really large x  <==>  2*pi*x  overflows
    let r = if Lrg_x {
        M_SQRT_2PI * x.sqrt() // sqrt(.): avoid overflow for very large x
    } else {
        M_2PI * x
    };
    if log {
        -yl - yh - (if Lrg_x { r.ln() } else { 0.5 * r.ln() })
    } else {
        (-yl).exp() * (-yh).exp() / (if Lrg_x { r } else { r.sqrt() })
    }
}

/// dpois() checks argument validity and calls dpois_raw().
pub fn dpois<R: Into<Real64>, P: Into<Positive64>>(x: R, lambda: P, log: bool) -> Real64 {
    let mut x = x.into().unwrap();
    let lambda = lambda.into().unwrap();

    if x.is_non_integer() {
        warn!("non-integer x = {}", x);
        return f64::d_0(log).into();
    }
    if x < 0.0 || !x.is_finite() {
        return f64::d_0(log).into();
    }

    x = x.round();

    dpois_raw(x, lambda, log).into()
}