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
108
109
110
111
112
113

// Translation of nmath's dbinom
//
// 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, PositiveInteger64, Probability64, Real64};

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

/// dbinom_raw() does the actual computation; note this is called by
/// other functions in addition to dbinom().
/// (1) dbinom_raw() has both p and q arguments, when one may be represented
/// more accurately than the other (in particular, in df()).
/// (2) dbinom_raw() does NOT check that inputs x and n are integers. This
/// should be done in the calling function, where necessary.
/// -- but is not the case at all when called e.g., from df() or dbeta() !
/// (3) Also does not check for 0 <= p <= 1 and 0 <= q <= 1 or NaN's.
/// Do this in the calling function.
pub fn dbinom_raw(x: f64, n: f64, p: f64, q: f64, log: bool) -> f64 {
    let mut lf = 0.0;
    let mut lc = 0.0;

    if p == 0.0 {
        return if x == 0.0 {
            f64::d_1(log)
        } else {
            f64::d_0(log)
        };
    }
    if q == 0.0 {
        return if x == n { f64::d_1(log) } else { f64::d_0(log) };
    }

    if x == 0.0 {
        if n == 0.0 {
            return f64::d_1(log);
        }
        lc = if p < 0.1 {
            -bd0(n, n * q).unwrap() - n * p
        } else {
            n * q.ln()
        };
        return lc.d_exp(log);
    }
    if x == n {
        lc = if q < 0.1 {
            -bd0(n, n * p).unwrap() - n * q
        } else {
            n * p.ln()
        };
        return lc.d_exp(log);
    }
    if x < 0.0 || x > n {
        return f64::d_0(log);
    }

    /* n*p or n*q can underflow to zero if n and p or q are small.  This
    used to occur in dbeta, and gives NaN as from R 2.3.0.0  */
    lc = stirlerr(n)
        - stirlerr(x)
        - stirlerr(n - x)
        - bd0(x, n * p).unwrap()
        - bd0(n - x, n * q).unwrap();

    /* f = (M_2PI*x*(n-x))/n; could overflow or underflow */
    /* Upto R 2.7.1:
     * lf = M_2PI.ln() + x.ln() + (n-x).ln() - n.ln();
     * -- following is much better for  x << n : */
    lf = strafe_consts::LN_2TPI + x.ln() + (-x / n).ln_1p();

    return (lc - 0.5 * lf).d_exp(log);
}

/// To compute the binomial probability, call dbinom(x,n,p).
/// This checks for argument validity, and calls dbinom_raw().
pub fn dbinom<R: Into<Real64>, PO: Into<PositiveInteger64>, PR: Into<Probability64>>(
    x: R,
    n: PO,
    p: PR,
    log: bool,
) -> Real64 {
    let mut x = x.into().unwrap();
    let mut n = n.into().unwrap();
    let p = p.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();
    }

    n = n.round();
    x = x.round();

    dbinom_raw(x, n, p, 1.0 - p, log).into()
}