Commit 0e3559e3 authored by Radford Neal's avatar Radford Neal

more work on Dpbinom

parent 434e0205
......@@ -146,7 +146,7 @@ cosh, sinh, tanh, acosh, asinh, atanh \cr
gamma, lgamma, digamma, trigamma, beta, lbeta \cr
\cr
dbeta \cr
dbinom \cr
dbinom, pbinom \cr
dcauchy, pcauchy, qcauchy, rcauchy \cr
dexp, pexp, qexp, rexp \cr
dgeom, pgeom \cr
......
......@@ -2930,7 +2930,7 @@ static void Ddbinom (double x, double n, double p,
if (dx || dn) abort();
if (!dp) return;
if (p <= 0 || p >= 1 || x < 0 || x > n) {
if (n <= 0 || p <= 0 || p >= 1 || x < 0 || x > n) {
*dp = 0;
}
else {
......@@ -2946,12 +2946,30 @@ static void Dpbinom (double q, double n, double p,
if (dq || dn) abort();
if (!dp) return;
if (p <= 0 || p >= 1 || q < 0 || q > n)
if (n <= 0 || p <= 0 || p >= 1 || q < 0 || q > n)
*dp = 0;
else {
*dp = n * dbinom(q,n-1,p,0);
if (lower_tail) *dp = - *dp;
if (log_p) *dp /= v;
/* Derivative wrt p of tail sum of binomial (n,p) probabilities gives
sum of pairs of binomial (n-1,p) probabilities, times +-n, with
all but one term then cancelling. Eg:
d/dp [ P(0;3,p) + P(1;3,p) + P(2;3,p) ] =
d/dp [ (1-p)^3 + 3*p*(1-p)^2 + 3*p^2*(1-p) ]
= 0 + 3*(1-p)^2 + 3*2*p*(1-p)
- 3*(1-p)^2 - 3*2*p*(1-p) - 3*p^2
= - 3*p^2
= - 3*P(2;2,p)
*/
if (log_p) {
double t = dbinom(q,n-1,p,TRUE) + log(n);
*dp = lower_tail ? -exp(t-v) : exp(t-v);
}
else {
double t = dbinom(q,n-1,p,FALSE)*n;
*dp = lower_tail ? -t : t;
}
}
}
......
......@@ -395,6 +395,11 @@ test3w(dbeta,log=TRUE)
test3v(dbinom)
test3v(dbinom,log=TRUE)
test3v(pbinom)
test3v(pbinom,log=TRUE)
test3v(pbinom,lower=FALSE)
test3v(pbinom,log=TRUE,lower=FALSE)
test3(dcauchy)
test3(dcauchy,log=TRUE)
test3z(dcauchy)
......
......@@ -921,6 +921,23 @@ r2 0.1429615 -0.9814315
r1 -1.94518 -6.865007
r2 -1.94518 -6.865006
>
> test3v(pbinom)
[,1] [,2]
r1 0.1977083 -1.639937
r2 0.1977083 -1.639937
> test3v(pbinom,log=TRUE)
[,1] [,2]
r1 -1.620963 -8.29473
r2 -1.620963 -8.29473
> test3v(pbinom,lower=FALSE)
[,1] [,2]
r1 0.8022917 1.639937
r2 0.8022917 1.639937
> test3v(pbinom,log=TRUE,lower=FALSE)
[,1] [,2]
r1 -0.220283 2.044065
r2 -0.220283 2.044065
>
> test3(dcauchy)
[,1] [,2]
r1 0.3414601 0.4543405
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment