Commit 41b3b34f authored by Radford Neal's avatar Radford Neal

implement derivatives for pcauchy and qcauchy

parent 8c0d67d3
......@@ -153,7 +153,7 @@ runif \cr
dexp, rexp \cr
dgeom \cr
dpois \cr
rcauchy \cr
pcauchy, qcauchy, rcauchy \cr
dnorm, pnorm, qnorm, rnorm \cr
rlnorm \cr
dlogis, plogis, qlogis, rlogis \cr
......
......@@ -2767,6 +2767,61 @@ SEXP do_Math2(SEXP call, SEXP op, SEXP args, SEXP env, int variant)
/* Derivatives of math3 functions. */
static void Dpcauchy (double q, double location, double scale,
double *dq, double *dlocation, double *dscale,
int lower_tail, int log_p, double v)
{
if (scale <= 0) {
if (dq) *dq = 0;
if (dlocation) *dlocation = 0;
if (dscale) *dscale = 0;
}
else {
double q0 = (q - location) / scale;
double dp0 = dcauchy (q0, 0, 1, 0);
if (dq) *dq = dp0 / scale;
if (dlocation) *dlocation = - dp0 / scale;
if (dscale) *dscale = - dp0 * q0 / scale;
if (!lower_tail) {
if (dq) *dq = -*dq;
if (dlocation) *dlocation = -*dlocation;
if (dscale) *dscale = -*dscale;
}
if (log_p) {
double expv = exp(-v);
if (dq) *dq *= expv;
if (dlocation) *dlocation *= expv;
if (dscale) *dscale *= expv;
}
}
}
static void Dqcauchy (double p, double location, double scale,
double *dp, double *dlocation, double *dscale,
int lower_tail, int log_p, double v)
{
if (scale <= 0) {
if (dp) *dp = 0;
if (dlocation) *dlocation = 0;
if (dscale) *dscale = 0;
}
else {
double q0 = (v - location) / scale;
if (dp) *dp = (lower_tail ? 1 : -1) * scale / dcauchy (q0, 0, 1, 0);
if (dlocation) *dlocation = 1;
if (dscale) *dscale = q0;
if (log_p)
if (dp) *dp *= exp(p);
}
}
static void Ddnorm (double x, double mu, double sigma,
double *dx, double *dmu, double *dsigma,
int give_log, double v)
......@@ -2937,8 +2992,8 @@ static struct { double (*fncall)(); void (*Dcall)(); } math3_table[48] = {
{ pbinom, 0 },
{ qbinom, 0 /* discrete */ },
{ dcauchy, 0 },
{ pcauchy, 0 },
{ qcauchy, 0 },
{ pcauchy, Dpcauchy },
{ qcauchy, Dqcauchy },
{ df, 0 },
{ pf, 0 },
{ qf, 0 },
......
......@@ -324,6 +324,20 @@ test2i(dgeom,log=TRUE)
test2i(dpois)
test2i(dpois,log=TRUE)
test3(pcauchy)
test3(pcauchy,log=TRUE)
test3(pcauchy,lower=FALSE)
test3(pcauchy,log=TRUE,lower=FALSE)
test3z(pcauchy)
test3z(pcauchy,log=TRUE)
test3z(pcauchy,lower=FALSE)
test3z(pcauchy,log=TRUE,lower=FALSE)
test3(qcauchy)
test3y(qcauchy,log=TRUE)
test3(qcauchy,lower=FALSE)
test3y(qcauchy,log=TRUE,lower=FALSE)
test2r(rcauchy)
test2r(rlnorm)
......
......@@ -652,6 +652,164 @@ r2 0.04879083 0.114805
r1 -3.020213 2.353004
r2 -3.020213 2.353004
>
> test3(pcauchy)
[,1] [,2]
r1 0.3232967 0.3414601
r2 0.3232967 0.3414601
[,1] [,2]
r1 0.3232967 -0.3414601
r2 0.3232967 -0.3414601
[,1] [,2]
r1 0.3232967 0.2117686
r2 0.3232967 0.2117687
x1 x2 x3
r1 0.3232967 0.3414601 -0.3414601 0.2117686
r2 0.3232967 0.3414601 -0.3414601 0.2117687
> test3(pcauchy,log=TRUE)
[,1] [,2]
r1 -1.129185 1.056182
r2 -1.129185 1.056182
[,1] [,2]
r1 -1.129185 -1.056182
r2 -1.129185 -1.056182
[,1] [,2]
r1 -1.129185 0.6550288
r2 -1.129185 0.6550288
x1 x2 x3
r1 -1.129185 1.056182 -1.056182 0.6550288
r2 -1.129185 1.056182 -1.056182 0.6550288
> test3(pcauchy,lower=FALSE)
[,1] [,2]
r1 0.6767033 -0.3414601
r2 0.6767033 -0.3414601
[,1] [,2]
r1 0.6767033 0.3414601
r2 0.6767033 0.3414601
[,1] [,2]
r1 0.6767033 -0.2117686
r2 0.6767033 -0.2117687
x1 x2 x3
r1 0.6767033 -0.3414601 0.3414601 -0.2117686
r2 0.6767033 -0.3414601 0.3414601 -0.2117687
> test3(pcauchy,log=TRUE,lower=FALSE)
[,1] [,2]
r1 -0.3905223 -0.5045935
r2 -0.3905223 -0.5045935
[,1] [,2]
r1 -0.3905223 0.5045935
r2 -0.3905223 0.5045935
[,1] [,2]
r1 -0.3905223 -0.3129416
r2 -0.3905223 -0.3129417
x1 x2 x3
r1 -0.3905223 -0.5045935 0.5045935 -0.3129416
r2 -0.3905223 -0.5045935 0.5045935 -0.3129417
> test3z(pcauchy)
[,1] [,2]
r1 0.4216113 0.04475888
r2 0.4216113 0.04475888
[,1] [,2]
r1 0.4216113 -0.04475888
r2 0.4216113 -0.04475888
[,1] [,2]
r1 0.4216113 0.01125093
r2 0.4216113 0.01125093
z1 z2 z3
r1 0.4216113 0.04475888 -0.04475888 0.01125093
r2 0.4216113 0.04475888 -0.04475888 0.01125093
> test3z(pcauchy,log=TRUE)
[,1] [,2]
r1 -0.8636715 0.1061615
r2 -0.8636715 0.1061615
[,1] [,2]
r1 -0.8636715 -0.1061615
r2 -0.8636715 -0.1061615
[,1] [,2]
r1 -0.8636715 0.02668555
r2 -0.8636715 0.02668555
z1 z2 z3
r1 -0.8636715 0.1061615 -0.1061615 0.02668555
r2 -0.8636715 0.1061615 -0.1061615 0.02668555
> test3z(pcauchy,lower=FALSE)
[,1] [,2]
r1 0.5783887 -0.04475888
r2 0.5783887 -0.04475888
[,1] [,2]
r1 0.5783887 0.04475888
r2 0.5783887 0.04475888
[,1] [,2]
r1 0.5783887 -0.01125093
r2 0.5783887 -0.01125093
z1 z2 z3
r1 0.5783887 -0.04475888 0.04475888 -0.01125093
r2 0.5783887 -0.04475888 0.04475888 -0.01125093
> test3z(pcauchy,log=TRUE,lower=FALSE)
[,1] [,2]
r1 -0.5475092 -0.07738546
r2 -0.5475092 -0.07738546
[,1] [,2]
r1 -0.5475092 0.07738546
r2 -0.5475092 0.07738546
[,1] [,2]
r1 -0.5475092 -0.01945219
r2 -0.5475092 -0.01945219
z1 z2 z3
r1 -0.5475092 -0.07738546 0.07738546 -0.01945219
r2 -0.5475092 -0.07738546 0.07738546 -0.01945219
>
> test3(qcauchy)
[,1] [,2]
r1 0.8463711 2.125985
r2 0.8463711 2.125985
[,1] [,2]
r1 0.8463711 1
r2 0.8463711 1
[,1] [,2]
r1 0.8463711 -0.07181421
r2 0.8463711 -0.07181422
x1 x2 x3
r1 0.8463711 2.125985 1 -0.07181421
r2 0.8463711 2.125985 1 -0.07181422
> test3y(qcauchy,log=TRUE)
[,1] [,2]
r1 0.01345346 3.85699
r2 0.01345346 3.85699
[,1] [,2]
r1 0.01345346 1
r2 0.01345346 1
[,1] [,2]
r1 0.01345346 0.676322
r2 0.01345346 0.676322
y1 y2 y3
r1 0.01345346 3.85699 1 0.676322
r2 0.01345346 3.85699 1 0.676322
> test3(qcauchy,lower=FALSE)
[,1] [,2]
r1 0.9430689 -2.125985
r2 0.9430689 -2.125985
[,1] [,2]
r1 0.9430689 1
r2 0.9430689 1
[,1] [,2]
r1 0.9430689 0.07181421
r2 0.9430689 0.07181422
x1 x2 x3
r1 0.9430689 -2.125985 1 0.07181421
r2 0.9430689 -2.125985 1 0.07181422
> test3y(qcauchy,log=TRUE,lower=FALSE)
[,1] [,2]
r1 -1.639653 -3.85699
r2 -1.639653 -3.85699
[,1] [,2]
r1 -1.639653 1
r2 -1.639653 1
[,1] [,2]
r1 -1.639653 -0.676322
r2 -1.639653 -0.676322
y1 y2 y3
r1 -1.639653 -3.85699 1 -0.676322
r2 -1.639653 -3.85699 1 -0.676322
>
> test2r(rcauchy)
x1 x2
r1 -1.585831 1 -2.305762
......
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