1 /*
2 * AUTHOR
3 * Catherine Loader, catherine@research.bell-labs.com.
4 * October 23, 2000.
5 *
6 * Merge in to R and further tweaks :
7 * Copyright (C) 2000, The R Core Team
8 * Copyright (C) 2008, The R Foundation
9 *
10 * This program is free software; you can redistribute it and/or modify
11 * it under the terms of the GNU General Public License as published by
12 * the Free Software Foundation; either version 2 of the License, or
13 * (at your option) any later version.
14 *
15 * This program is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 * GNU General Public License for more details.
19 *
20 * You should have received a copy of the GNU General Public License
21 * along with this program; if not, a copy is available at
22 * http://www.r-project.org/Licenses/
23 *
24 *
25 * DESCRIPTION
26 *
27 * To compute the binomial probability, call dbinom(x,n,p).
28 * This checks for argument validity, and calls dbinom_raw().
29 *
30 * dbinom_raw() does the actual computation; note this is called by
31 * other functions in addition to dbinom().
32 * (1) dbinom_raw() has both p and q arguments, when one may be represented
33 * more accurately than the other (in particular, in df()).
34 * (2) dbinom_raw() does NOT check that inputs x and n are integers. This
35 * should be done in the calling function, where necessary.
36 * -- but is not the case at all when called e.g., from df() or dbeta() !
37 * (3) Also does not check for 0 <= p <= 1 and 0 <= q <= 1 or NaN's.
38 * Do this in the calling function.
39 */
40
41 #include "nmath.h"
42 #include "dpq.h"
43
44 double attribute_hidden
dbinom_raw(double x,double n,double p,double q,int give_log)45 dbinom_raw(double x, double n, double p, double q, int give_log)
46 {
47 double lf, lc;
48
49 if (p == 0) return((x == 0) ? R_D__1 : R_D__0);
50 if (q == 0) return((x == n) ? R_D__1 : R_D__0);
51
52 if (x == 0) {
53 if(n == 0) return R_D__1;
54 lc = (p < 0.1) ? -bd0(n,n*q) - n*p : n*log(q);
55 return( R_D_exp(lc) );
56 }
57 if (x == n) {
58 lc = (q < 0.1) ? -bd0(n,n*p) - n*q : n*log(p);
59 return( R_D_exp(lc) );
60 }
61 if (x < 0 || x > n) return( R_D__0 );
62
63 /* n*p or n*q can underflow to zero if n and p or q are small. This
64 used to occur in dbeta, and gives NaN as from R 2.3.0. */
65 lc = stirlerr(n) - stirlerr(x) - stirlerr(n-x) - bd0(x,n*p) - bd0(n-x,n*q);
66
67 /* f = (M_2PI*x*(n-x))/n; could overflow or underflow */
68 /* Upto R 2.7.1:
69 * lf = log(M_2PI) + log(x) + log(n-x) - log(n);
70 * -- following is much better for x << n : */
71 lf = log(M_2PI) + log(x) + log1p(- x/n);
72
73 return R_D_exp(lc - 0.5*lf);
74 }
75
dbinom(double x,double n,double p,int give_log)76 double dbinom(double x, double n, double p, int give_log)
77 {
78 #ifdef IEEE_754
79 /* NaNs propagated correctly */
80 if (ISNAN(x) || ISNAN(n) || ISNAN(p)) return x + n + p;
81 #endif
82
83 if (p < 0 || p > 1 || R_D_negInonint(n))
84 ML_ERR_return_NAN;
85 R_D_nonint_check(x);
86 if (x < 0 || !R_FINITE(x)) return R_D__0;
87
88 n = R_D_forceint(n);
89 x = R_D_forceint(x);
90
91 return dbinom_raw(x, n, p, 1-p, give_log);
92 }
93