Commit 59cf0f6e authored by Radford Neal's avatar Radford Neal

implement derivative for pgeom

parent e5741192
......@@ -151,7 +151,7 @@ gamma, lgamma, digamma, trigamma, beta, lbeta \cr
\cr
runif \cr
dexp, pexp, qexp, rexp \cr
dgeom \cr
dgeom, pgeom \cr
dpois \cr
pcauchy, qcauchy, rcauchy \cr
dnorm, pnorm, qnorm, rnorm \cr
......
......@@ -2573,7 +2573,7 @@ static void Ddgeom (double x, double p, double *dx /*ignored*/,
{
if (!dp) return;
if (x < 0 || p <= 0 || p >= 1)
if (x < 0 || p <= 0 || x > 0 && p >= 1)
*dp = 0;
else {
*dp = 1/p - x/(1-p);
......@@ -2582,6 +2582,24 @@ static void Ddgeom (double x, double p, double *dx /*ignored*/,
}
}
static void Dpgeom (double q, double p, double *dq /*ignored*/, double *dp,
int lower_tail, int log_p, double v)
{
if (!dp) return;
if (q < 0 || p <= 0)
*dp = 0;
else {
double t = (q+1) / (p-1);
if (log_p)
*dp = lower_tail ? -t*expm1(-v) : t;
else
*dp = lower_tail ? t*(v-1) : t*v;
}
}
static void Ddpois (double x, double lambda, double *dx /*ignored*/,
double *dlambda, double v, int give_log)
{
......@@ -2628,7 +2646,7 @@ static struct { double (*fncall)(); void (*Dcall)(); } math2_table[31] = {
{ pexp, Dpexp },
{ qexp, Dqexp },
{ dgeom, Ddgeom },
{ pgeom, 0 },
{ pgeom, Dpgeom },
{ qgeom, 0 /* discrete */ },
{ dpois, Ddpois },
{ ppois, 0 },
......
......@@ -348,6 +348,11 @@ test1r(rexp)
test2i(dgeom)
test2i(dgeom,log=TRUE)
test2i(pgeom)
test2i(pgeom,log=TRUE)
test2i(pgeom,lower=FALSE)
test2i(pgeom,log=TRUE,lower=FALSE)
test2i(dpois)
test2i(dpois,log=TRUE)
......
......@@ -766,6 +766,23 @@ r2 0.001044058 -0.02858399
r1 -6.86464 -27.37777
r2 -6.86464 -27.37777
>
> test2i(pgeom)
[,1] [,2]
r1 0.9998771 0.004667647
r2 0.9998771 0.004667643
> test2i(pgeom,log=TRUE)
[,1] [,2]
r1 -0.0001228599 0.004668215
r2 -0.0001228599 0.004668216
> test2i(pgeom,lower=FALSE)
[,1] [,2]
r1 0.0001228524 -0.004667642
r2 0.0001228524 -0.004667643
> test2i(pgeom,log=TRUE,lower=FALSE)
[,1] [,2]
r1 -9.004527 -37.99392
r2 -9.004527 -37.99392
>
> test2i(dpois)
[,1] [,2]
r1 0.04879083 0.114805
......
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