Commit 26282881 authored by Radford Neal's avatar Radford Neal

implement gradients for runif

parent d05e0046
......@@ -44,11 +44,11 @@ R_NORETURN static void invalid(SEXP call)
/* Random sampling from 1-parameter families. */
/* Derivatives for random1 generators. */
/* Derivatives for random1 generators of continuous distributions. */
static double Drexp (double r, double scale)
{
return r / scale;
return scale <= 0 ? 0 : r / scale;
}
/* Table of functions for generation and derivatives of random1 distibutions. */
......@@ -59,10 +59,10 @@ static struct {
{
{ rchisq, 0 },
{ rexp, Drexp },
{ rgeom, 0 },
{ rpois, 0 },
{ rgeom, 0 /* discrete */ },
{ rpois, 0 /* discrete */ },
{ rt, 0 },
{ rsignrank,0 }
{ rsignrank,0 /* discrete */ }
};
static SEXP do_random1(SEXP call, SEXP op, SEXP args, SEXP rho)
......@@ -172,13 +172,21 @@ static SEXP do_random1(SEXP call, SEXP op, SEXP args, SEXP rho)
/* Random sampling from 2-parameter families. */
/* Derivatives for random2 generators. */
/* Derivatives for random2 generators of continuous distributions. */
static void Drnorm (double r, double mu, double sigma,
double *dmu, double *dsigma)
{
if (dmu) *dmu = 1;
if (dsigma) *dsigma = sigma == 0 ? 0 : (r - mu) / sigma;
if (dsigma) *dsigma = sigma <= 0 ? 0 : (r - mu) / sigma;
}
static void Drunif (double r, double a, double b,
double *da, double *db)
{
double u = b <= a ? 0.5 : (r - a) / (b - a);
if (da) *da = 1-u;
if (db) *db = u;
}
static struct {
......@@ -187,19 +195,19 @@ static struct {
} rand2_table[14] =
{
{ rbeta, 0 },
{ rbinom, 0 },
{ rbinom, 0 /* discrete */ },
{ rcauchy, 0 },
{ rf, 0 },
{ rgamma, 0 },
{ rlnorm, 0 },
{ rlogis, 0 },
{ rnbinom, 0 },
{ rnbinom, 0 /* discrete */ },
{ rnorm, Drnorm },
{ runif, 0 },
{ runif, Drunif },
{ rweibull, 0 },
{ rwilcox, 0 },
{ rwilcox, 0 /* discrete */ },
{ rnchisq, 0 },
{ rnbinom_mu, 0 }
{ rnbinom_mu, 0 /* discrete */ }
};
static SEXP do_random2(SEXP call, SEXP op, SEXP args, SEXP rho)
......@@ -326,14 +334,17 @@ static SEXP do_random2(SEXP call, SEXP op, SEXP args, SEXP rho)
}
/* Random sampling from 3-parameter families. */
/* Random sampling from 3-parameter families. There's only one, and it's
a discrete distribution, so there are currently no derivatives to compute,
and hence the code to do so has not been written. (It would follow the
pattern in do_random2, with one more parameter.) */
static struct {
double (*fncall)(double, double, double);
void (*Dcall)(double, double, double, double, double *, double *, double *);
} rand3_table[1] =
{
{ rhyper, 0 }
{ rhyper, 0 /* discrete */ }
};
static SEXP do_random3(SEXP call, SEXP op, SEXP args, SEXP rho)
......@@ -948,24 +959,24 @@ attribute_hidden FUNTAB R_FunTab_random[] =
{"rchisq", do_random1, 0, 51000011, 2, {PP_FUNCALL, PREC_FN, 0}},
{"rexp", do_random1, 1, 51000011, 2, {PP_FUNCALL, PREC_FN, 0}},
{"rgeom", do_random1, 2, 51000011, 2, {PP_FUNCALL, PREC_FN, 0}},
{"rpois", do_random1, 3, 51000011, 2, {PP_FUNCALL, PREC_FN, 0}},
{"rgeom", do_random1, 2, 1000011, 2, {PP_FUNCALL, PREC_FN, 0}},
{"rpois", do_random1, 3, 1000011, 2, {PP_FUNCALL, PREC_FN, 0}},
{"rt", do_random1, 4, 51000011, 2, {PP_FUNCALL, PREC_FN, 0}},
{"rsignrank", do_random1, 5, 1000011, 2, {PP_FUNCALL, PREC_FN, 0}},
{"rbeta", do_random2, 0, 1000011, 3, {PP_FUNCALL, PREC_FN, 0}},
{"rbeta", do_random2, 0, 61000011, 3, {PP_FUNCALL, PREC_FN, 0}},
{"rbinom", do_random2, 1, 1000011, 3, {PP_FUNCALL, PREC_FN, 0}},
{"rcauchy", do_random2, 2, 1000011, 3, {PP_FUNCALL, PREC_FN, 0}},
{"rf", do_random2, 3, 1000011, 3, {PP_FUNCALL, PREC_FN, 0}},
{"rgamma", do_random2, 4, 1000011, 3, {PP_FUNCALL, PREC_FN, 0}},
{"rlnorm", do_random2, 5, 1000011, 3, {PP_FUNCALL, PREC_FN, 0}},
{"rlogis", do_random2, 6, 1000011, 3, {PP_FUNCALL, PREC_FN, 0}},
{"rcauchy", do_random2, 2, 61000011, 3, {PP_FUNCALL, PREC_FN, 0}},
{"rf", do_random2, 3, 61000011, 3, {PP_FUNCALL, PREC_FN, 0}},
{"rgamma", do_random2, 4, 61000011, 3, {PP_FUNCALL, PREC_FN, 0}},
{"rlnorm", do_random2, 5, 61000011, 3, {PP_FUNCALL, PREC_FN, 0}},
{"rlogis", do_random2, 6, 61000011, 3, {PP_FUNCALL, PREC_FN, 0}},
{"rnbinom", do_random2, 7, 1000011, 3, {PP_FUNCALL, PREC_FN, 0}},
{"rnbinom_mu", do_random2, 13, 1000011, 3, {PP_FUNCALL, PREC_FN, 0}},
{"rnchisq", do_random2, 12, 1000011, 3, {PP_FUNCALL, PREC_FN, 0}},
{"rnorm", do_random2, 8, 61000011, 3, {PP_FUNCALL, PREC_FN, 0}},
{"runif", do_random2, 9, 1000011, 3, {PP_FUNCALL, PREC_FN, 0}},
{"rweibull", do_random2, 10, 1000011, 3, {PP_FUNCALL, PREC_FN, 0}},
{"runif", do_random2, 9, 61000011, 3, {PP_FUNCALL, PREC_FN, 0}},
{"rweibull", do_random2, 10, 61000011, 3, {PP_FUNCALL, PREC_FN, 0}},
{"rwilcox", do_random2, 11, 1000011, 3, {PP_FUNCALL, PREC_FN, 0}},
{"rnchisq", do_random2, 12, 61000011, 3, {PP_FUNCALL, PREC_FN, 0}},
{"rnbinom_mu", do_random2, 13, 1000011, 3, {PP_FUNCALL, PREC_FN, 0}},
{"rhyper", do_random3, 0, 1000011, 4, {PP_FUNCALL, PREC_FN, 0}},
{"sample", do_sample, 0, 1000011, 4, {PP_FUNCALL, PREC_FN, 0}},
{"rmultinom", do_rmultinom, 0, 1000011, 3, {PP_FUNCALL, PREC_FN, 0}},
......
......@@ -140,7 +140,7 @@ print (with_gradient (a) sin(a))
# Check consistency of results between with_gradient and numericDeriv.
x <- 0.32739
x1 <- 0.89472; x2 <- 0.49718
x1 <- 0.47718; x2 <- 0.89472
i1 <- 3
bindgrads <- function (r1,r2)
......@@ -229,3 +229,5 @@ test2i(dgeom)
test2i(dgeom,log=TRUE)
test2r(rnorm)
test2r(runif)
......@@ -249,7 +249,7 @@ attr(,"gradient")
> # Check consistency of results between with_gradient and numericDeriv.
>
> x <- 0.32739
> x1 <- 0.89472; x2 <- 0.49718
> x1 <- 0.47718; x2 <- 0.89472
> i1 <- 3
>
> bindgrads <- function (r1,r2)
......@@ -402,44 +402,44 @@ r2 10.43204 -58.12982
>
> test2(atan2)
[,1] [,2]
r1 1.063601 0.4745389
r2 1.063601 0.4745389
r1 0.4899538 0.8701601
r2 0.4899538 0.8701601
[,1] [,2]
r1 1.063601 -0.8539753
r2 1.063601 -0.8539753
r1 0.4899538 -0.4640815
r2 0.4899538 -0.4640815
x1 x2
r1 1.063601 0.4745389 -0.8539753
r2 1.063601 0.4745389 -0.8539753
r1 0.4899538 0.8701601 -0.4640815
r2 0.4899538 0.8701601 -0.4640815
g1 g2
r1 1.063601 0.4745389 -0.8539753
r2 1.063601 0.4745389 -0.8539753
r1 0.4899538 0.8701601 -0.4640815
r2 0.4899538 0.8701601 -0.4640815
>
> test2(dexp)
[,1] [,2]
r1 0.318657 -0.1584299
r2 0.318657 -0.1584299
r1 0.5838058 -0.5223427
r2 0.5838058 -0.5223427
[,1] [,2]
r1 0.318657 0.35582
r2 0.318657 0.35582
r1 0.5838058 0.3739206
r2 0.5838058 0.3739206
x1 x2
r1 0.318657 -0.1584299 0.35582
r2 0.318657 -0.1584299 0.35582
r1 0.5838058 -0.5223427 0.3739206
r2 0.5838058 -0.5223427 0.3739206
g1 g2
r1 0.318657 -0.1584299 0.35582
r2 0.318657 -0.1584299 0.35582
r1 0.5838058 -0.5223427 0.3739206
r2 0.5838058 -0.5223427 0.3739206
> test2(dexp,log=TRUE)
[,1] [,2]
r1 -1.14364 -0.49718
r2 -1.14364 -0.49718
r1 -0.5381869 -0.89472
r2 -0.5381869 -0.89472
[,1] [,2]
r1 -1.14364 1.116624
r2 -1.14364 1.116624
r1 -0.5381869 0.6404881
r2 -0.5381869 0.6404881
x1 x2
r1 -1.14364 -0.49718 1.116624
r2 -1.14364 -0.49718 1.116624
r1 -0.5381869 -0.89472 0.6404881
r2 -0.5381869 -0.89472 0.6404881
g1 g2
r1 -1.14364 -0.49718 1.116624
r2 -1.14364 -0.49718 1.116624
r1 -0.5381869 -0.89472 0.6404881
r2 -0.5381869 -0.89472 0.6404881
> test1r(rexp)
[,1] [,2]
r1 0.7957275 -2.430519
......@@ -447,15 +447,20 @@ r2 0.7957275 -2.430519
>
> test2i(dgeom)
[,1] [,2]
r1 0.06320498 -0.2499760
r2 0.06320498 -0.2499761
r1 0.001044058 -0.02858399
r2 0.001044058 -0.02858399
> test2i(dgeom,log=TRUE)
[,1] [,2]
r1 -2.761372 -3.955006
r2 -2.761372 -3.955006
r1 -6.86464 -27.37777
r2 -6.86464 -27.37777
>
> test2r(rnorm)
x1 x2
r1 1.060049 1 0.3325331
r2 1.060049 1 0.3325331
r1 0.774704 1 0.3325331
r2 0.774704 1 0.3325331
>
> test2r(runif)
x1 x2
r1 0.7403373 0.3697434 0.6302566
r2 0.7403373 0.3697434 0.6302566
>
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