[Desire-commits] r9 - in packages: . loglognorm loglognorm/R loglognorm/man loglognorm/src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu May 15 00:55:15 CEST 2008
Author: olafm
Date: 2008-05-15 00:55:15 +0200 (Thu, 15 May 2008)
New Revision: 9
Added:
packages/loglognorm/
packages/loglognorm/DESCRIPTION
packages/loglognorm/NAMESPACE
packages/loglognorm/R/
packages/loglognorm/R/loglognorm.R
packages/loglognorm/man/
packages/loglognorm/man/dloglognorm.Rd
packages/loglognorm/src/
packages/loglognorm/src/loglognorm.c
Log:
Inital checkin of loglognorm package
Added: packages/loglognorm/DESCRIPTION
===================================================================
--- packages/loglognorm/DESCRIPTION (rev 0)
+++ packages/loglognorm/DESCRIPTION 2008-05-14 22:55:15 UTC (rev 9)
@@ -0,0 +1,14 @@
+Package: loglognorm
+Version: 0.9.0
+Date: 2008-05-10
+Title: Double log normal distribution functions
+Author:
+ Heike Trautmann <trautmann at statistik.uni-dortmund.de>
+ and Detlef Steuer <detlef.steuer at hsu-hamburg.de>
+ and Olaf Mersmann <olafm at statistik.uni-dortmund.de>
+Maintainer: desiRe developers <desire-developers at lists.r-forge.r-project.org>
+Depends: R (>= 2.7.0)
+Description:
+License: GPLv2
+URL: http://r-forge.r-project.org/projects/desire
+LazyData: yes
Added: packages/loglognorm/NAMESPACE
===================================================================
--- packages/loglognorm/NAMESPACE (rev 0)
+++ packages/loglognorm/NAMESPACE 2008-05-14 22:55:15 UTC (rev 9)
@@ -0,0 +1,7 @@
+useDynLib(loglognorm)
+
+export(dloglognorm)
+export(ploglognorm)
+export(qloglognorm)
+export(rloglognorm)
+export(eloglognorm)
Added: packages/loglognorm/R/loglognorm.R
===================================================================
--- packages/loglognorm/R/loglognorm.R (rev 0)
+++ packages/loglognorm/R/loglognorm.R 2008-05-14 22:55:15 UTC (rev 9)
@@ -0,0 +1,25 @@
+##
+## loglognorm.R - Interface to loglognorm.c
+##
+## Authors:
+## Heike Trautmann <trautmann at statistik.uni-dortmund.de>
+## Detlef Steuer <detlef.steuer at hsu-hamburg.de>
+## Olaf Mersmann <olafm at statistik.uni-dortmund.de>
+##
+
+dloglognorm <- function(x, mean=0, sd=1)
+ .External("dloglognorm", x, mean, sd, PACKAGE="loglognorm")
+
+ploglognorm <- function(q, mean=0, sd=1)
+ .External("ploglognorm", q, mean, sd, PACKAGE="loglognorm")
+
+qloglognorm <- function(p, mean=0, sd=1)
+ .External("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")
+
+eloglognorm <- function(mean=0, sd=1)
+ .External("eloglognorm", mean, sd, PACKAGE="loglognorm")
+
Added: packages/loglognorm/man/dloglognorm.Rd
===================================================================
--- packages/loglognorm/man/dloglognorm.Rd (rev 0)
+++ packages/loglognorm/man/dloglognorm.Rd 2008-05-14 22:55:15 UTC (rev 9)
@@ -0,0 +1,43 @@
+\name{dloglognorm}
+\alias{dloglognorm}
+\alias{ploglognorm}
+\alias{qloglognorm}
+\alias{rloglognorm}
+\alias{eloglognorm}
+\title{The Double Log Normal Distribution}
+\description{
+ Density, distribution function, quantile function, random generation
+ and expected value function for the double log normal distribution
+ with mean equal to 'mean' and standard deviation equal to 'sd'.
+}
+\usage{
+dloglognorm(x, mean = 0, sd = 1)
+ploglognorm(q, mean = 0, sd = 1)
+qloglognorm(p, mean = 0, sd = 1)
+rloglognorm(n, mean = 0, sd = 1)
+eloglognorm(mean = 0, sd = 1)
+}
+\arguments{
+ \item{x,q}{vector of quantiles.}
+ \item{p}{vector of probabilites.}
+ \item{n}{number of observations.}
+ \item{mean}{vector of means.}
+ \item{sd}{vector of standard deviations.}
+}
+\details{
+ If 'mean' or 'sd' are not specified they assume the default values of
+ '0' and '1', respectively.
+}
+\value{
+ 'dloglognorm' gives the density, 'ploglognorm' gives the distribution
+ function, 'qloglognorm' gives the quantile function, 'rloglognorm'
+ generates random deviates and 'eloglognorm' gives the expected value
+ of the distirbution.
+}
+\references{ FIXME: Erstes auftauchen? }
+\author{
+ Heike Trautmann \email{trautmann at statistik.uni-dortmund.de},
+ Detlef Steuer \email{steuer at hsu-hamburg.de} and
+ Olaf Mersmann \email{olafm at statistik.uni-dortmund.de}
+}
+\keyword{distribution}
Added: packages/loglognorm/src/loglognorm.c
===================================================================
--- packages/loglognorm/src/loglognorm.c (rev 0)
+++ packages/loglognorm/src/loglognorm.c 2008-05-14 22:55:15 UTC (rev 9)
@@ -0,0 +1,183 @@
+/*
+ * loglognorm.c - Implementation of double log normal functions
+ *
+ * Authors:
+ * Heike Trautmann <trautmann at statistik.uni-dortmund.de>
+ * Detlef Steuer <detlef.steuer at hsu-hamburg.de>
+ * Olaf Mersmann <olafm at statistik.uni-dortmund.de>
+ */
+
+#include <R.h>
+#include <Rmath.h>
+#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);
+
+#define UNPACK_REAL_VALUE(ARGS, SXP, DBL) \
+ ARGS = CDR(ARGS); \
+ SXP = CAR(ARGS); \
+ PROTECT(SXP = coerceVector(SXP, REALSXP)); \
+ DBL = *REAL(SXP);
+
+#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) \
+ PROTECT(SXP = allocVector(REALSXP, SIZE)); \
+ DBL = REAL(SXP);
+
+
+SEXP dloglognorm(SEXP args) {
+ 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(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);
+
+ 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);
+ ret[i] = c3 * exp(c4);
+ } else {
+ ret[i] = 0.0;
+ }
+ }
+ UNPROTECT(4);
+ return s_ret;
+}
+
+
+SEXP ploglognorm(SEXP args) {
+ R_len_t i, n;
+ SEXP s_q, s_mean, s_sd, s_ret;
+ double *q, mean, sd, *ret;
+
+ 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);
+
+ for (i = 0; i < n; ++i) {
+ if (q[i] < 0.0) {
+ ret[i] = 0.0;
+ } else if (q[i] > 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);
+ }
+ }
+ UNPROTECT(4);
+ return s_ret;
+}
+
+
+SEXP qloglognorm(SEXP args) {
+ R_len_t i, n;
+ SEXP s_p, s_mean, s_sd, s_ret;
+ double *p, mean, sd, *ret;
+
+ UNPACK_REAL_VECTOR(args, s_p, p);
+ UNPACK_REAL_VALUE(args, s_mean, mean);
+ UNPACK_REAL_VALUE(args, s_sd, sd);
+
+ 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;
+ }
+ }
+ UNPROTECT(4);
+ return s_ret;
+}
+
+
+typedef struct {
+ double mean, sd;
+} loglognorm_param;
+
+static void loglognorm_intgr(double *x, int n, void *ex) {
+ int i;
+ loglognorm_param *lp = (loglognorm_param *)ex;
+ const double mean = lp->mean;
+ const double sd = lp->sd;
+
+ /* Taken from Trautmann (2004) p. 54 */
+ for (i = 0; i < n; ++i) {
+ x[i] = exp(-exp(mean + sd * x[i])) * M_1_SQRT_2PI * exp(-0.5 * pow(x[i], 2.0));
+ }
+}
+
+SEXP eloglognorm(SEXP args) {
+ R_len_t i, n_mean, n_sd;
+ SEXP s_mean, s_sd, s_ret;
+ double *mean, *sd, *ret, tmp;
+
+ UNPACK_REAL_VECTOR(args, s_mean, mean);
+ UNPACK_REAL_VECTOR(args, s_sd, sd);
+
+ n_mean = length(s_mean);
+ n_sd = length(s_sd);
+
+ if (n_mean != n_sd)
+ error("Length of mean and sd differ (%i != %i)", n_mean, n_sd);
+
+ ALLOC_REAL_VECTOR(n_mean, s_ret, ret);
+
+ /* 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;
+ 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];
+ Rdqagi(loglognorm_intgr, (void *)&lp, &bound, &inf,
+ &epsabs, &epsrel,
+ &result, &abserr, &neval, &ier,
+ &limit, &lenw, &last, iwork, work);
+ /* FIXME: Possibly check agains lower bound given in Trautmann
+ * (2004):
+ *
+ * E(X) >= exp(-1.0) * Phi(-mean/sd)
+ */
+ if (ier >= 1) { /* Failure */
+ ret[i] = R_NaN;
+ } else { /* No error */
+ ret[i] = result;
+ }
+ }
+ UNPROTECT(3);
+ return s_ret;
+}
+
More information about the Desire-commits
mailing list