[Desire-commits] r20 - in packages/loglognorm: . R man src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed May 21 14:10:50 CEST 2008
Author: olafm
Date: 2008-05-21 14:10:49 +0200 (Wed, 21 May 2008)
New Revision: 20
Modified:
packages/loglognorm/DESCRIPTION
packages/loglognorm/R/loglognorm.R
packages/loglognorm/man/dloglognorm.Rd
packages/loglognorm/src/loglognorm.c
Log:
* Fixed usage of pnorm(), dnorm() and qnorm()
* Reworked functions so that they can deal with arguments of arbitrary (different) lenghts.
Modified: packages/loglognorm/DESCRIPTION
===================================================================
--- packages/loglognorm/DESCRIPTION 2008-05-21 09:24:28 UTC (rev 19)
+++ packages/loglognorm/DESCRIPTION 2008-05-21 12:10:49 UTC (rev 20)
@@ -1,6 +1,6 @@
Package: loglognorm
-Version: 0.9.5
-Date: 2008-05-19
+Version: 0.9.6
+Date: 2008-05-20
Title: Double log normal distribution functions
Author:
Heike Trautmann <trautmann at statistik.uni-dortmund.de>
Modified: packages/loglognorm/R/loglognorm.R
===================================================================
--- packages/loglognorm/R/loglognorm.R 2008-05-21 09:24:28 UTC (rev 19)
+++ packages/loglognorm/R/loglognorm.R 2008-05-21 12:10:49 UTC (rev 20)
@@ -8,23 +8,23 @@
##
dloglognorm <- function(x, mean=0, sd=1)
- .External("dloglognorm", x, mean, sd, PACKAGE="loglognorm")
+ .Call("dloglognorm", x, mean, sd, PACKAGE="loglognorm")
ploglognorm <- function(q, mean=0, sd=1)
- .External("ploglognorm", q, mean, sd, PACKAGE="loglognorm")
+ .Call("ploglognorm", q, mean, sd, PACKAGE="loglognorm")
qloglognorm <- function(p, mean=0, sd=1)
- .External("qloglognorm", p, mean, sd, PACKAGE="loglognorm")
+ .Call("qloglognorm", p, mean, sd, PACKAGE="loglognorm")
## FIXME: Faster way?
rloglognorm <- function(n, mean=0, sd=1)
- .External("qloglognorm", runif(n), mean=mean, sd=sd, PACKAGE="loglognorm")
+ .Call("qloglognorm", runif(n), mean=mean, sd=sd, PACKAGE="loglognorm")
mloglognorm <- function(moment, mean, sd)
- .External("mloglognorm", mean, sd, moment, PACKAGE="loglognorm")
+ .Call("mloglognorm", mean, sd, moment, PACKAGE="loglognorm")
eloglognorm <- function(mean, sd)
- .External("mloglognorm", mean, sd, rep(1, length(mean)), PACKAGE="loglognorm")
+ .Call("mloglognorm", mean, sd, rep(1, length(mean)), PACKAGE="loglognorm")
vloglognorm <- function(mean, sd) {
m1 <- mloglognorm(rep(1, length(mean)), mean, sd)
Modified: packages/loglognorm/man/dloglognorm.Rd
===================================================================
--- packages/loglognorm/man/dloglognorm.Rd 2008-05-21 09:24:28 UTC (rev 19)
+++ packages/loglognorm/man/dloglognorm.Rd 2008-05-21 12:10:49 UTC (rev 20)
@@ -49,4 +49,16 @@
Detlef Steuer \email{steuer at hsu-hamburg.de} and
Olaf Mersmann \email{olafm at statistik.uni-dortmund.de}
}
+\examples{
+ x <- seq(0, 1, by=0.05)
+ ## Several different shapes of the density:
+ par(mfrow=c(3, 1))
+ curve(dloglognorm(x, -0.2, 0.2), 0, 1, main="DLN(-0.2, 0.2)")
+ curve(dloglognorm(x, 0.2, 1.0), 0, 1, main="DLN(0.2, 2.0)")
+ curve(dloglognorm(x, 0.2, 1.8), 0, 1, main="DLN(0.2, 2.0)")
+
+ ## Check precision:
+ z <- x - pnorm(qnorm(x, .2, 1.0), .2, 1.0)
+ max(z)
+}
\keyword{distribution}
Modified: packages/loglognorm/src/loglognorm.c
===================================================================
--- packages/loglognorm/src/loglognorm.c 2008-05-21 09:24:28 UTC (rev 19)
+++ packages/loglognorm/src/loglognorm.c 2008-05-21 12:10:49 UTC (rev 20)
@@ -12,104 +12,96 @@
#include <Rinternals.h>
#include <R_ext/Applic.h>
-#define UNPACK_REAL_VECTOR(ARGS, SXP, DBL) \
- ARGS = CDR(ARGS); \
- SXP = CAR(ARGS); \
- PROTECT(SXP = coerceVector(SXP, REALSXP)); \
- DBL = REAL(SXP);
+#ifndef MAX
+#define MAX(A, B) ((A > B)?(A):(B))
+#endif
-#define UNPACK_REAL_VALUE(ARGS, SXP, DBL) \
- ARGS = CDR(ARGS); \
- SXP = CAR(ARGS); \
- PROTECT(SXP = coerceVector(SXP, REALSXP)); \
- DBL = *REAL(SXP);
+#define UNPACK_REAL_VECTOR(S, D, N) \
+ double *D = REAL(S); \
+ R_len_t N = length(S);
-#define UNPACK_INTEGER_VALUE(ARGS, SXP, INT) \
- ARGS = CDR(ARGS); \
- SXP = CAR(ARGS); \
- PROTECT(SXP = coerceVector(SXP, INTSXP)); \
- INT = *INTEGER(SXP);
-
-#define ALLOC_REAL_VECTOR(SIZE, SXP, DBL) \
+#define ALLOC_REAL_VECTOR(SXP, DBL, SIZE) \
+ SEXP SXP; \
PROTECT(SXP = allocVector(REALSXP, SIZE)); \
- DBL = REAL(SXP);
+ double *DBL = REAL(SXP);
-SEXP dloglognorm(SEXP args) {
+SEXP dloglognorm(SEXP s_x, SEXP s_mean, SEXP s_sd) {
R_len_t i, n;
- SEXP s_x, s_mean, s_sd, s_ret;
- double *x, mean, sd, *ret, c1, c2, c3, c4;
+ UNPACK_REAL_VECTOR(s_x , x , n_x);
+ UNPACK_REAL_VECTOR(s_mean, mean, n_mean);
+ UNPACK_REAL_VECTOR(s_sd , sd , n_sd);
- UNPACK_REAL_VECTOR(args, s_x, x);
- UNPACK_REAL_VALUE(args, s_mean, mean);
- UNPACK_REAL_VALUE(args, s_sd, sd);
-
- n = length(s_x);
- ALLOC_REAL_VECTOR(n, s_ret, ret);
+ n = MAX(MAX(n_x, n_mean), n_sd);
+ ALLOC_REAL_VECTOR(s_ret, ret, n);
- c1 = - M_1_SQRT_2PI / sd;
- c2 = -1.0/(2.0 * sd * sd);
for (i = 0; i < n; ++i) {
- if (0 <= x[i] && x[i] <= 1) {
- c3 = c1 / (log(x[i]) * x[i]);
- c4 = c2 * pow(log(-log(x[i])) - mean, 2);
+ const double cx = x[i % n_x];
+ if (0 <= cx && cx <= 1) {
+ const double cmean = mean[i % n_mean];
+ const double csd = sd[i % n_sd];
+
+ const double c1 = -1.0/(2.0 * csd * csd);
+ const double c2 = log(cx);
+ const double c3 = -M_1_SQRT_2PI / (csd * c2 * cx);
+ const double c4 = c1 * pow(log(-c2) - cmean, 2);
ret[i] = c3 * exp(c4);
} else {
ret[i] = 0.0;
}
}
- UNPROTECT(4);
+ UNPROTECT(1); /* s_ret */
return s_ret;
}
-SEXP ploglognorm(SEXP args) {
+SEXP ploglognorm(SEXP s_q, SEXP s_mean, SEXP s_sd) {
R_len_t i, n;
- SEXP s_q, s_mean, s_sd, s_ret;
- double *q, mean, sd, *ret;
+ UNPACK_REAL_VECTOR(s_q , q , n_q);
+ UNPACK_REAL_VECTOR(s_mean, mean, n_mean);
+ UNPACK_REAL_VECTOR(s_sd , sd , n_sd);
- UNPACK_REAL_VECTOR(args, s_q, q);
- UNPACK_REAL_VALUE(args, s_mean, mean);
- UNPACK_REAL_VALUE(args, s_sd, sd);
-
- n = length(s_q);
- ALLOC_REAL_VECTOR(n, s_ret, ret);
+ n = MAX(MAX(n_q, n_mean), n_sd);
+ ALLOC_REAL_VECTOR(s_ret, ret, n);
for (i = 0; i < n; ++i) {
- if (q[i] < 0.0) {
+ const double cq = q[i % n_q];
+ const double cmean = mean[i % n_mean];
+ const double csd = sd[i % n_sd];
+ if (cq < 0.0) {
ret[i] = 0.0;
- } else if (q[i] > 1.0) {
+ } else if (cq > 1.0) {
ret[i] = 1.0;
- } else {
- /* FIXME: Better use log=TRUE? */
- ret[i] = 1.0 - pnorm((log(-log(q[i])) - mean)/sd, 0.0, 1.0, FALSE, FALSE);
+ } else { /* q \in [0, 1] */
+ /* Directly return the upper tail: */
+ ret[i] = pnorm(log(-log(cq)), cmean, csd, FALSE, FALSE);
}
}
- UNPROTECT(4);
+ UNPROTECT(1); /* s_ret */
return s_ret;
}
-SEXP qloglognorm(SEXP args) {
+SEXP qloglognorm(SEXP s_p, SEXP s_mean, SEXP s_sd) {
R_len_t i, n;
- SEXP s_p, s_mean, s_sd, s_ret;
- double *p, mean, sd, *ret;
+ UNPACK_REAL_VECTOR(s_p , p , n_p);
+ UNPACK_REAL_VECTOR(s_mean, mean, n_mean);
+ UNPACK_REAL_VECTOR(s_sd , sd , n_sd);
- UNPACK_REAL_VECTOR(args, s_p, p);
- UNPACK_REAL_VALUE(args, s_mean, mean);
- UNPACK_REAL_VALUE(args, s_sd, sd);
+ n = MAX(MAX(n_p, n_mean), n_sd);
+ ALLOC_REAL_VECTOR(s_ret, ret, n);
- n = length(s_p);
- ALLOC_REAL_VECTOR(n, s_ret, ret);
-
for (i = 0; i < n; ++i) {
- if (0 <= p[i] && p[i] <= 1.0) {
- ret[i] = exp(-exp(sd * qnorm(1-p[i], 0.0, 1.0, FALSE, FALSE) + mean));
- } else {
- ret[i] = R_NaN;
+ const double cp = p[i % n_p];
+ const double cmean = mean[i % n_mean];
+ const double csd = sd[i % n_sd];
+ if (0 <= cp && cp <= 1.0) { /* p \in [0, 1] */
+ ret[i] = exp(-exp(qnorm(1-cp, cmean, csd, TRUE, FALSE)));
+ } else {
+ ret[i] = 0.0;
}
}
- UNPROTECT(4);
+ UNPROTECT(1); /* s_ret */
return s_ret;
}
@@ -131,39 +123,32 @@
}
}
-SEXP mloglognorm(SEXP args) {
- R_len_t i, n_mean, n_sd, n_moment;
- SEXP s_mean, s_sd, s_moment, s_ret;
- double *mean, *sd, *moment, *ret, tmp;
+SEXP mloglognorm(SEXP s_moment, SEXP s_mean, SEXP s_sd) {
+ R_len_t i, n;
+ UNPACK_REAL_VECTOR(s_moment, moment, n_moment);
+ UNPACK_REAL_VECTOR(s_mean , mean , n_mean);
+ UNPACK_REAL_VECTOR(s_sd , sd , n_sd);
- UNPACK_REAL_VECTOR(args, s_mean, mean);
- UNPACK_REAL_VECTOR(args, s_sd, sd);
- UNPACK_REAL_VECTOR(args, s_moment, moment);
- n_mean = length(s_mean);
- n_sd = length(s_sd);
- n_moment = length(s_moment);
- if (n_mean != n_sd)
- error("Length of mean, sd and moment differ (%i != %i)", n_mean, n_sd);
-
- ALLOC_REAL_VECTOR(n_mean, s_ret, ret);
-
- /* Parameters for Rdqagi() */
+ n = MAX(MAX(n_moment, n_mean), n_sd);
+ ALLOC_REAL_VECTOR(s_ret, ret, n);
+
+ /* Parameters for Rdqagi: */
int limit = 100;
int lenw = 4 * limit;
int *iwork = (int *) R_alloc(limit, sizeof(int));
double *work = (double *) R_alloc(lenw, sizeof(double));
loglognorm_param lp;
double bound = 0; /* not used */
- int inf = 2; /* 2 = -1nf - inf */
- double epsabs = 0.0000001;
+ int inf = 2; /* 2 = -Inf - inf */
+ double epsabs = 0.00000001;
double epsrel = epsabs;
double result, abserr;
int neval, ier, last;
for (i = 0; i < n_mean; ++i) {
- lp.mean = mean[i];
- lp.sd = sd[i];
- lp.r = moment[i];
+ lp.mean = mean[i % n_mean];
+ lp.sd = sd[i % n_sd];
+ lp.r = moment[i % n_moment];
Rdqagi(loglognorm_intgr, (void *)&lp, &bound, &inf,
&epsabs, &epsrel,
&result, &abserr, &neval, &ier,
@@ -179,7 +164,7 @@
ret[i] = result;
}
}
- UNPROTECT(4); /* s_mean, s_sd, s_moment, s_ret */
+ UNPROTECT(1); /* s_ret */
return s_ret;
}
More information about the Desire-commits
mailing list