Commit e5741192 authored by Radford Neal's avatar Radford Neal

implement derivatives for pexp and qexp

parent 0d9e4d40
......@@ -150,7 +150,7 @@ cosh, sinh, tanh, acosh, asinh, atanh \cr
gamma, lgamma, digamma, trigamma, beta, lbeta \cr
\cr
runif \cr
dexp, rexp \cr
dexp, pexp, qexp, rexp \cr
dgeom \cr
dpois \cr
pcauchy, qcauchy, rcauchy \cr
......
......@@ -2459,6 +2459,9 @@ SEXP do_abs(SEXP call, SEXP op, SEXP args, SEXP env, int variant)
}
/* MATHEMATICAL FUNCTIONS OF TWO NUMERIC ARGUMENTS
(plus 0, 1, or 2 integers). */
/* Derivatives of math2 functions. */
static void Datan2 (double y, double x, double *dy, double *dx, double v)
......@@ -2518,6 +2521,53 @@ static void Ddexp (double x, double scale, double *dx, double *dscale,
}
}
static void Dpexp (double q, double scale, double *dq, double *dscale,
int lower_tail, int log_p, double v)
{
if (scale <= 0) {
if (dq) *dq = 0;
if (dscale) *dscale = 0;
}
else {
double q0 = q / scale;
double dp0 = dexp (q0, 1, 0);
if (dq) *dq = dp0 / scale;
if (dscale) *dscale = - dp0 * q0 / scale;
if (!lower_tail) {
if (dq) *dq = -*dq;
if (dscale) *dscale = -*dscale;
}
if (log_p) {
double expv = exp(-v);
if (dq) *dq *= expv;
if (dscale) *dscale *= expv;
}
}
}
static void Dqexp (double p, double scale, double *dp, double *dscale,
int lower_tail, int log_p, double v)
{
if (scale <= 0) {
if (dp) *dp = 0;
if (dscale) *dscale = 0;
}
else {
double q0 = v / scale;
if (dp) *dp = (lower_tail ? 1 : -1) * scale / dexp (q0, 1, 0);
if (dscale) *dscale = q0;
if (log_p)
if (dp) *dp *= exp(p);
}
}
static void Ddgeom (double x, double p, double *dx /*ignored*/,
double *dp, double v, int give_log)
{
......@@ -2575,8 +2625,8 @@ static struct { double (*fncall)(); void (*Dcall)(); } math2_table[31] = {
{ pchisq, 0 },
{ qchisq, 0 },
{ dexp, Ddexp },
{ pexp, 0 },
{ qexp, 0 },
{ pexp, Dpexp },
{ qexp, Dqexp },
{ dgeom, Ddgeom },
{ pgeom, 0 },
{ qgeom, 0 /* discrete */ },
......@@ -2598,8 +2648,6 @@ static struct { double (*fncall)(); void (*Dcall)(); } math2_table[31] = {
{ fprec, 0 }
};
/* Mathematical functions of 2 numeric arguments (plus 0, 1, or 2 integers) */
SEXP do_math2 (SEXP call, SEXP op, SEXP args, SEXP env)
{
SEXP res;
......@@ -2765,6 +2813,9 @@ SEXP do_Math2(SEXP call, SEXP op, SEXP args, SEXP env, int variant)
}
/* MATHEMATICAL FUNCTIONS OF THREE NUMERIC ARGUMENTS
(plus 0, 1, or 2 integers). */
/* Derivatives of math3 functions. */
static void Dpcauchy (double q, double location, double scale,
......@@ -3034,8 +3085,6 @@ static struct { double (*fncall)(); void (*Dcall)(); } math3_table[48] = {
{ qnbinom_mu, 0 /* discrete */ }
};
/* Mathematical functions of 3 numeric arguments (plus 0, 1, or 2 integers) */
SEXP do_math3 (SEXP call, SEXP op, SEXP args, SEXP env)
{
SEXP res;
......@@ -3160,7 +3209,7 @@ SEXP do_math3 (SEXP call, SEXP op, SEXP args, SEXP env)
}
/* Mathematical Functions of Four (Real) Arguments */
/* MATHEMATICAL FUNCTIONS OF FOUR (REAL) ARGUMENTS. */
static void setup_Math4 (SEXP *sa, SEXP *sb, SEXP *sc, SEXP *sd, SEXP *sy,
int na, int nb, int nc, int nd, SEXP lcall)
......
......@@ -212,6 +212,22 @@ test2 <- function (fun,...) {
))
}
test2y <- function (fun,...) {
print (bindgrads (numericDeriv(quote(fun(y1,y3,...)),"y1"),
with_gradient (y1) fun(y1,y3,...)))
print (bindgrads (numericDeriv(quote(fun(y1,y3,...)),"y3"),
with_gradient (y3) fun(y1,y3,...)))
print (bindgrads (numericDeriv(quote(fun(y1,y3,...)),c("y1","y3")),
with_gradient (y1,y3) fun(y1,y3,...)))
print (bindgrads (numericDeriv(quote(fun(y1,y3,...)),c("y1","y3")),
{ r <- with_gradient (y1) { s <- with_gradient (y3) fun(y1,y3,...);
g2 <<- attr(s,"gradient"); s }
attr(r,"gradient") <- cbind(g1=attr(r,"gradient"),g2=g2)
r
}
))
}
test2z <- function (fun,...) {
print (bindgrads (numericDeriv(quote(fun(z1,z2,...)),"z1"),
with_gradient (z1) fun(z1,z2,...)))
......@@ -316,6 +332,17 @@ test2z(lbeta)
test2(dexp)
test2(dexp,log=TRUE)
test2(pexp)
test2(pexp,log=TRUE)
test2(pexp,lower=FALSE)
test2(pexp,log=TRUE,lower=FALSE)
test2(qexp)
test2y(qexp,log=TRUE)
test2(qexp,lower=FALSE)
test2y(qexp,log=TRUE,lower=FALSE)
test1r(rexp)
test2i(dgeom)
......
......@@ -357,6 +357,22 @@ attr(,"gradient")
+ ))
+ }
>
> test2y <- function (fun,...) {
+ print (bindgrads (numericDeriv(quote(fun(y1,y3,...)),"y1"),
+ with_gradient (y1) fun(y1,y3,...)))
+ print (bindgrads (numericDeriv(quote(fun(y1,y3,...)),"y3"),
+ with_gradient (y3) fun(y1,y3,...)))
+ print (bindgrads (numericDeriv(quote(fun(y1,y3,...)),c("y1","y3")),
+ with_gradient (y1,y3) fun(y1,y3,...)))
+ print (bindgrads (numericDeriv(quote(fun(y1,y3,...)),c("y1","y3")),
+ { r <- with_gradient (y1) { s <- with_gradient (y3) fun(y1,y3,...);
+ g2 <<- attr(s,"gradient"); s }
+ attr(r,"gradient") <- cbind(g1=attr(r,"gradient"),g2=g2)
+ r
+ }
+ ))
+ }
>
> test2z <- function (fun,...) {
+ print (bindgrads (numericDeriv(quote(fun(z1,z2,...)),"z1"),
+ with_gradient (z1) fun(z1,z2,...)))
......@@ -629,6 +645,113 @@ r2 -0.5381869 -0.89472 0.6404881
g1 g2
r1 -0.5381869 -0.89472 0.6404881
r2 -0.5381869 -0.89472 0.6404881
>
> test2(pexp)
[,1] [,2]
r1 0.3474989 0.5838058
r2 0.3474989 0.5838058
[,1] [,2]
r1 0.3474989 0.3113605
r2 0.3474989 0.3113605
x1 x2
r1 0.3474989 0.5838058 0.3113605
r2 0.3474989 0.5838058 0.3113605
g1 g2
r1 0.3474989 0.5838058 0.3113605
r2 0.3474989 0.5838058 0.3113605
> test2(pexp,log=TRUE)
[,1] [,2]
r1 -1.056994 1.680022
r2 -1.056994 1.680022
[,1] [,2]
r1 -1.056994 0.8960041
r2 -1.056994 0.8960041
x1 x2
r1 -1.056994 1.680022 0.8960041
r2 -1.056994 1.680022 0.8960041
g1 g2
r1 -1.056994 1.680022 0.8960041
r2 -1.056994 1.680022 0.8960041
> test2(pexp,lower=FALSE)
[,1] [,2]
r1 0.6525011 -0.5838057
r2 0.6525011 -0.5838058
[,1] [,2]
r1 0.6525011 -0.3113605
r2 0.6525011 -0.3113605
x1 x2
r1 0.6525011 -0.5838057 -0.3113605
r2 0.6525011 -0.5838058 -0.3113605
g1 g2
r1 0.6525011 -0.5838057 -0.3113605
r2 0.6525011 -0.5838058 -0.3113605
> test2(pexp,log=TRUE,lower=FALSE)
[,1] [,2]
r1 -0.4269425 -0.89472
r2 -0.4269425 -0.89472
[,1] [,2]
r1 -0.4269425 -0.47718
r2 -0.4269425 -0.47718
x1 x2
r1 -0.4269425 -0.89472 -0.47718
r2 -0.4269425 -0.89472 -0.47718
g1 g2
r1 -0.4269425 -0.89472 -0.47718
r2 -0.4269425 -0.89472 -0.47718
>
> test2(qexp)
[,1] [,2]
r1 0.7248279 2.137768
r2 0.7248279 2.137768
[,1] [,2]
r1 0.7248279 -0.810117
r2 0.7248279 -0.810117
x1 x2
r1 0.7248279 2.137768 -0.810117
r2 0.7248279 2.137768 -0.810117
g1 g2
r1 0.7248279 2.137768 -0.810117
r2 0.7248279 2.137768 -0.810117
> test2y(qexp,log=TRUE)
[,1] [,2]
r1 0.9564287 1.815181
r2 0.9564287 1.815181
[,1] [,2]
r1 0.9564287 -0.7825916
r2 0.9564287 -0.7825916
y1 y3
r1 0.9564287 1.815181 -0.7825916
r2 0.9564287 1.815181 -0.7825916
g1 g2
r1 0.9564287 1.815181 -0.7825916
r2 0.9564287 1.815181 -0.7825916
> test2(qexp,lower=FALSE)
[,1] [,2]
r1 0.8269196 -2.342236
r2 0.8269196 -2.342236
[,1] [,2]
r1 0.8269196 -0.9242216
r2 0.8269196 -0.9242217
x1 x2
r1 0.8269196 -2.342236 -0.9242216
r2 0.8269196 -2.342236 -0.9242217
g1 g2
r1 0.8269196 -2.342236 -0.9242216
r2 0.8269196 -2.342236 -0.9242217
> test2y(qexp,log=TRUE,lower=FALSE)
[,1] [,2]
r1 0.3044684 -0.8182436
r2 0.3044684 -0.8182436
[,1] [,2]
r1 0.3044684 -0.2491293
r2 0.3044684 -0.2491293
y1 y3
r1 0.3044684 -0.8182436 -0.2491293
r2 0.3044684 -0.8182436 -0.2491293
g1 g2
r1 0.3044684 -0.8182436 -0.2491293
r2 0.3044684 -0.8182436 -0.2491293
>
> test1r(rexp)
[,1] [,2]
r1 0.7957275 -2.430519
......
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