| Line 4... |
Line 4... |
| 4 |
* October 23, 2000 and Feb, 2001.
|
4 |
* October 23, 2000 and Feb, 2001.
|
| 5 |
*
|
5 |
*
|
| 6 |
* dnbinom_mu(): Martin Maechler, June 2008
|
6 |
* dnbinom_mu(): Martin Maechler, June 2008
|
| 7 |
*
|
7 |
*
|
| 8 |
* Merge in to R and improvements notably for |x| << size :
|
8 |
* Merge in to R and improvements notably for |x| << size :
|
| 9 |
* Copyright (C) 2000--2021, The R Core Team
|
9 |
* Copyright (C) 2000--2025, The R Core Team
|
| 10 |
*
|
10 |
*
|
| 11 |
* This program is free software; you can redistribute it and/or modify
|
11 |
* This program is free software; you can redistribute it and/or modify
|
| 12 |
* it under the terms of the GNU General Public License as published by
|
12 |
* it under the terms of the GNU General Public License as published by
|
| 13 |
* the Free Software Foundation; either version 2 of the License, or
|
13 |
* the Free Software Foundation; either version 2 of the License, or
|
| 14 |
* (at your option) any later version.
|
14 |
* (at your option) any later version.
|
| Line 53... |
Line 53... |
| 53 |
return(give_log ? size*log(prob) : pow(prob, size));
|
53 |
return(give_log ? size*log(prob) : pow(prob, size));
|
| 54 |
}
|
54 |
}
|
| 55 |
if(!R_FINITE(size)) size = DBL_MAX;
|
55 |
if(!R_FINITE(size)) size = DBL_MAX;
|
| 56 |
|
56 |
|
| 57 |
if(x < 1e-10 * size) { // instead of dbinom_raw(), use 2 terms of Abramowitz & Stegun (6.1.47)
|
57 |
if(x < 1e-10 * size) { // instead of dbinom_raw(), use 2 terms of Abramowitz & Stegun (6.1.47)
|
| - |
|
58 |
double xx2s = /* x(x-1)/(2*size) robustly */
|
| - |
|
59 |
(x < sqrt(DBL_MAX)) ? ldexp(x*(x-1), -1)/size : x*(ldexp(x,-1)/size);
|
| 58 |
return R_D_exp(size * log(prob) + x * (log(size) + log1p(-prob))
|
60 |
return R_D_exp(size * log(prob) + x * (log(size) + log1p(-prob))
|
| 59 |
- lgamma1p(x) + log1p(x*(x-1)/(2*size)));
|
61 |
- lgamma1p(x) + log1p(xx2s));
|
| 60 |
} else {
|
62 |
} else {
|
| 61 |
/* log( size/(size+x) ) is much less accurate than log1p(- x/(size+x))
|
63 |
/* log( size/(size+x) ) is much less accurate than log1p(- x/(size+x))
|
| 62 |
for |x| << size (and actually when x < size): */
|
64 |
for |x| << size (and actually when x < size): */
|
| 63 |
double p = give_log ? (x < size ? log1p(-x/(size+x)) : log(size/(size+x)))
|
65 |
double p = give_log ? (x < size ? log1p(-x/(size+x)) : log(size/(size+x)))
|
| 64 |
: size/(size+x),
|
66 |
: size/(size+x),
|
| Line 94... |
Line 96... |
| 94 |
if(x == 0)/* be accurate, both for n << mu, and n >> mu :*/
|
96 |
if(x == 0)/* be accurate, both for n << mu, and n >> mu :*/
|
| 95 |
return R_D_exp(size * (size < mu ? log(size/(size+mu)) : log1p(- mu/(size+mu))));
|
97 |
return R_D_exp(size * (size < mu ? log(size/(size+mu)) : log1p(- mu/(size+mu))));
|
| 96 |
if(x < 1e-10 * size) { /* don't use dbinom_raw() but MM's formula: */
|
98 |
if(x < 1e-10 * size) { /* don't use dbinom_raw() but MM's formula: */
|
| 97 |
/* FIXME --- 1e-8 shows problem; rather use algdiv() from ./toms708.c */
|
99 |
/* FIXME --- 1e-8 shows problem; rather use algdiv() from ./toms708.c */
|
| 98 |
double p = (size < mu ? log(size/(1 + size/mu)) : log(mu / (1 + mu/size)));
|
100 |
double p = (size < mu ? log(size/(1 + size/mu)) : log(mu / (1 + mu/size)));
|
| 99 |
return R_D_exp(x * p - mu - lgamma1p(x) +
|
101 |
double xx2s = /* x(x-1)/(2*size) robustly */
|
| - |
|
102 |
(x < sqrt(DBL_MAX)) ? ldexp(x*(x-1), -1)/size : x*(ldexp(x,-1)/size);
|
| 100 |
log1p(x*(x-1)/(2*size)));
|
103 |
return R_D_exp(x * p - mu - lgamma1p(x) + log1p(xx2s));
|
| 101 |
} else {
|
104 |
} else {
|
| 102 |
/* no unnecessary cancellation inside dbinom_raw, when
|
105 |
/* no unnecessary cancellation inside dbinom_raw, when
|
| 103 |
x_ = size and n_ = x+size are so close that n_ - x_ loses accuracy
|
106 |
x_ = size and n_ = x+size are so close that n_ - x_ loses accuracy
|
| 104 |
but log( size/(size+x) ) is much less accurate than log1p(- x/(size+x))
|
107 |
but log( size/(size+x) ) is much less accurate than log1p(- x/(size+x))
|
| 105 |
for |x| << size (and actually when x < size): */
|
108 |
for |x| << size (and actually when x < size): */
|
The R Project SVN
R