/src/r-source/src/nmath/dpois.c
Line | Count | Source |
1 | | /* |
2 | | * AUTHOR |
3 | | * Catherine Loader, catherine@research.bell-labs.com. |
4 | | * October 23, 2000. |
5 | | * |
6 | | * Merge in to R: |
7 | | * Copyright (C) 2000-2021 The R Core Team |
8 | | * |
9 | | * This program is free software; you can redistribute it and/or modify |
10 | | * it under the terms of the GNU General Public License as published by |
11 | | * the Free Software Foundation; either version 2 of the License, or |
12 | | * (at your option) any later version. |
13 | | * |
14 | | * This program is distributed in the hope that it will be useful, |
15 | | * but WITHOUT ANY WARRANTY; without even the implied warranty of |
16 | | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
17 | | * GNU General Public License for more details. |
18 | | * |
19 | | * You should have received a copy of the GNU General Public License |
20 | | * along with this program; if not, a copy is available at |
21 | | * https://www.R-project.org/Licenses/ |
22 | | * |
23 | | * |
24 | | * DESCRIPTION |
25 | | * |
26 | | * dpois() checks argument validity and calls dpois_raw(). |
27 | | * |
28 | | * dpois_raw() computes the Poisson probability lb^x exp(-lb) / x!. |
29 | | * This does not check that x is an integer, since dgamma() may |
30 | | * call this with a fractional x argument. Any necessary argument |
31 | | * checks should be done in the calling function. |
32 | | * |
33 | | */ |
34 | | |
35 | | #include "nmath.h" |
36 | | #include "dpq.h" |
37 | | |
38 | 0 | #define M_SQRT_2PI 2.50662827463100050241576528481104525301 /* sqrt(2*pi) */ |
39 | | // sqrt(2 * Rmpfr::Const("pi", 128)) |
40 | 0 | #define x_LRG 2.86111748575702815380240589208115399625e+307 /* = 2^1023 / pi */ |
41 | | |
42 | | // called also from dgamma.c, pgamma.c, dnbeta.c, dnbinom.c, dnchisq.c : |
43 | | double dpois_raw(double x, double lambda, int give_log) |
44 | 0 | { |
45 | | /* x >= 0 ; integer for dpois(), but not e.g. for pgamma()! |
46 | | lambda >= 0 |
47 | | */ |
48 | 0 | if (lambda == 0) return( (x == 0) ? R_D__1 : R_D__0 ); |
49 | 0 | if (!R_FINITE(lambda)) return R_D__0; // including for the case where x = lambda = +Inf |
50 | 0 | if (x < 0) return( R_D__0 ); |
51 | 0 | if (x <= lambda * DBL_MIN) return(R_D_exp(-lambda) ); |
52 | 0 | if (lambda < x * DBL_MIN) { |
53 | 0 | if (!R_FINITE(x)) // lambda < x = +Inf |
54 | 0 | return R_D__0; |
55 | | // else |
56 | 0 | return(R_D_exp(-lambda + x*log(lambda) -lgammafn(x+1))); |
57 | 0 | } |
58 | | // R <= 4.0.x had return(R_D_fexp( M_2PI*x, -stirlerr(x)-bd0(x,lambda) )); |
59 | 0 | double yh, yl; |
60 | 0 | ebd0 (x, lambda, &yh, &yl); |
61 | 0 | yl += stirlerr(x); |
62 | 0 | bool Lrg_x = (x >= x_LRG); //really large x <==> 2*pi*x overflows |
63 | 0 | double r = Lrg_x |
64 | 0 | ? M_SQRT_2PI * sqrt(x) // sqrt(.): avoid overflow for very large x |
65 | 0 | : M_2PI * x; |
66 | 0 | return give_log |
67 | 0 | ? -yl - yh - (Lrg_x ? log(r) : 0.5 * log(r)) |
68 | 0 | : exp(-yl) * exp(-yh) / (Lrg_x ? r : sqrt(r)); |
69 | 0 | } |
70 | | |
71 | | double dpois(double x, double lambda, int give_log) |
72 | 0 | { |
73 | 0 | #ifdef IEEE_754 |
74 | 0 | if(ISNAN(x) || ISNAN(lambda)) |
75 | 0 | return x + lambda; |
76 | 0 | #endif |
77 | | |
78 | 0 | if (lambda < 0) ML_WARN_return_NAN; |
79 | 0 | R_D_nonint_check(x); |
80 | 0 | if (x < 0 || !R_FINITE(x)) |
81 | 0 | return R_D__0; |
82 | | |
83 | 0 | x = R_forceint(x); |
84 | |
|
85 | 0 | return( dpois_raw(x,lambda,give_log) ); |
86 | 0 | } |