[Desire-commits] r10 - in packages: . truncnorm truncnorm/R truncnorm/man truncnorm/src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu May 15 00:55:57 CEST 2008
Author: olafm
Date: 2008-05-15 00:55:57 +0200 (Thu, 15 May 2008)
New Revision: 10
Added:
packages/truncnorm/
packages/truncnorm/DESCRIPTION
packages/truncnorm/NAMESPACE
packages/truncnorm/R/
packages/truncnorm/R/truncnorm.R
packages/truncnorm/man/
packages/truncnorm/man/dtruncnorm.Rd
packages/truncnorm/src/
packages/truncnorm/src/truncnorm.c
Log:
Initial checkin of truncnorm package
Added: packages/truncnorm/DESCRIPTION
===================================================================
--- packages/truncnorm/DESCRIPTION (rev 0)
+++ packages/truncnorm/DESCRIPTION 2008-05-14 22:55:57 UTC (rev 10)
@@ -0,0 +1,14 @@
+Package: truncnorm
+Version: 0.9.0
+Date: 2008-05-12
+Title: Trunctated Normal Distribution
+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/truncnorm/NAMESPACE
===================================================================
--- packages/truncnorm/NAMESPACE (rev 0)
+++ packages/truncnorm/NAMESPACE 2008-05-14 22:55:57 UTC (rev 10)
@@ -0,0 +1,7 @@
+useDynLib(truncnorm)
+
+export(dtruncnorm)
+export(ptruncnorm)
+export(qtruncnorm)
+export(rtruncnorm)
+export(etruncnorm)
Added: packages/truncnorm/R/truncnorm.R
===================================================================
--- packages/truncnorm/R/truncnorm.R (rev 0)
+++ packages/truncnorm/R/truncnorm.R 2008-05-14 22:55:57 UTC (rev 10)
@@ -0,0 +1,46 @@
+##
+## truncnorm.R - Interface to truncnorm.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>
+##
+
+dtruncnorm <- function(x, a, b, mean=0, sd=1)
+ .External("dtruncnorm", x, a, b, mean, sd)
+
+ptruncnorm <- function(q, a, b, mean=0, sd=1)
+ .External("ptruncnorm", q, a, b, mean, sd)
+
+qtruncnorm <- function(p, a, b, mean=0, sd=1) {
+ ## Taken out of package 'msm'
+ ret <- numeric(length(p))
+ ret[p == 1] <- Inf
+ ret[p == 0] <- -Inf
+ ret[p < 0 | p > 1] <- NaN
+ ind <- (p > 0 & p < 1)
+ if (any(ind)) {
+ hind <- seq(along = p)[ind]
+ h <- function(y) {
+ (ptruncnorm(y, a, b, mean, sd) - p)[hind[i]]
+ }
+ ptmp <- numeric(length(p[ind]))
+ for (i in 1:length(p[ind])) {
+ interval <- c(-1, 1)
+ while (h(interval[1]) * h(interval[2]) >= 0) interval <- interval +
+ c(-1, 1) * 0.5 * (interval[2] - interval[1])
+ ptmp[i] <- uniroot(h, interval, tol = .Machine$double.eps)$root
+ }
+ ret[ind] <- ptmp
+ }
+ return (ret)
+}
+
+rtruncnorm <- function(n, a, b, mean=0, sd=1)
+ .External("rtruncnorm", n, a, b, mean, sd)
+
+etruncnorm <- function(a, b, mean=0, sd=1)
+ .External("etruncnorm", a, b, mean, sd)
+
+
Added: packages/truncnorm/man/dtruncnorm.Rd
===================================================================
--- packages/truncnorm/man/dtruncnorm.Rd (rev 0)
+++ packages/truncnorm/man/dtruncnorm.Rd 2008-05-14 22:55:57 UTC (rev 10)
@@ -0,0 +1,47 @@
+\name{dtruncnorm}
+\alias{dtruncnorm}
+\alias{ptruncnorm}
+\alias{qtruncnorm}
+\alias{rtruncnorm}
+\alias{etruncnorm}
+\title{The Truncated Normal Distribution}
+\description{
+ Density, distribution function, quantile function, random generation
+ and expected value function for the truncated normal distribution
+ with mean equal to 'mean' and standard deviation equal to 'sd'.
+}
+\usage{
+dtruncnorm(x, a, b, mean = 0, sd = 1)
+ptruncnorm(q, a, b, mean = 0, sd = 1)
+qtruncnorm(p, a, b, mean = 0, sd = 1)
+rtruncnorm(n, a, b, mean = 0, sd = 1)
+etruncnorm(a, b, mean = 0, sd = 1)
+}
+\arguments{
+ \item{x,q}{vector of quantiles.}
+ \item{p}{vector of probabilites.}
+ \item{n}{number of observations.}
+ \item{a}{vector of lower bounds.}
+ \item{b}{vector of upper bounds.}
+ \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{
+ 'dtruncnorm' gives the density, 'ptruncnorm' gives the distribution
+ function, 'qtruncnorm' gives the quantile function, 'rtruncnorm'
+ generates random deviates and 'etruncnorm' 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}.
+ 'qtruncnorm' was originally published in package 'msm' written by
+ C. H. Jackson \email{chris.jackson at mrc-bsu.cam.ac.uk}
+}
+\keyword{distribution}
Added: packages/truncnorm/src/truncnorm.c
===================================================================
--- packages/truncnorm/src/truncnorm.c (rev 0)
+++ packages/truncnorm/src/truncnorm.c 2008-05-14 22:55:57 UTC (rev 10)
@@ -0,0 +1,155 @@
+/*
+ * truncnorm.c - Implementation of truncated normal distribution
+ *
+ * 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 dtruncnorm(SEXP args) {
+ R_len_t i, n;
+ SEXP s_x, s_a, s_b, s_mean, s_sd, s_ret;
+ double *x, a, b, mean, sd, c1, c2, *ret;
+
+ UNPACK_REAL_VECTOR(args, s_x, x);
+ UNPACK_REAL_VALUE(args, s_a, a);
+ UNPACK_REAL_VALUE(args, s_b, b);
+ 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);
+
+ c2 = sd * pnorm(b, mean, sd, FALSE, FALSE) - pnorm(a, mean, sd, FALSE, FALSE);
+ for (i = 0; i < n; ++i) {
+ double t = x[i];
+ if (a <= t && t <= b) { /* In range: */
+ c1 = dnorm(t, mean, sd, FALSE);
+ ret[i] = - c1 / c2;
+ } else { /* Truncated: */
+ ret[i] = 0.0;
+ }
+ }
+ UNPROTECT(6);
+ return s_ret;
+}
+
+
+SEXP ptruncnorm(SEXP args) {
+ R_len_t i, n;
+ SEXP s_q, s_a, s_b, s_mean, s_sd, s_ret;
+ double *q, a, b, mean, sd, cx, ca, cb, *ret;
+
+ UNPACK_REAL_VECTOR(args, s_q, q);
+ UNPACK_REAL_VALUE(args, s_a, a);
+ UNPACK_REAL_VALUE(args, s_b, b);
+ 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);
+
+ ca = pnorm(a, mean, sd, FALSE, FALSE);
+ cb = pnorm(b, mean, sd, FALSE, FALSE);
+
+ for (i = 0; i < n; ++i) {
+ if (q[i] < a) {
+ ret[i] = 0.0;
+ } else if (q[i] > b) {
+ ret[i] = 1.0;
+ } else {
+ cx = pnorm(q[i], mean, sd, FALSE, FALSE);
+ ret[i] = (cx - ca) / (cb - ca);
+ }
+ }
+ UNPROTECT(6);
+ return s_ret;
+}
+
+
+SEXP rtruncnorm(SEXP args) {
+ R_len_t i;
+ int n;
+ SEXP s_n, s_a, s_b, s_mean, s_sd, s_ret;
+ double a, b, mean, sd, *ret, tmp;
+
+ UNPACK_INTEGER_VALUE(args, s_n, n);
+ UNPACK_REAL_VALUE(args, s_a, a);
+ UNPACK_REAL_VALUE(args, s_b, b);
+ UNPACK_REAL_VALUE(args, s_mean, mean);
+ UNPACK_REAL_VALUE(args, s_sd, sd);
+
+ ALLOC_REAL_VECTOR(n, s_ret, ret);
+ GetRNGstate();
+ for (i = 0; i < n; ++i) {
+ while (1) {
+ tmp = rnorm(mean, sd);
+ if (a <= tmp && tmp <= b)
+ break;
+ }
+ ret[i] = tmp;
+ }
+ PutRNGstate();
+ UNPROTECT(6);
+ return s_ret;
+}
+
+
+SEXP etruncnorm(SEXP args) {
+ R_len_t i, na, nb;
+ SEXP s_a, s_b, s_mean, s_sd, s_ret;
+ double *a, *b, mean, sd, *ret;
+
+ UNPACK_REAL_VECTOR(args, s_a, a);
+ UNPACK_REAL_VECTOR(args, s_b, b);
+ UNPACK_REAL_VALUE(args, s_mean, mean);
+ UNPACK_REAL_VALUE(args, s_sd, sd);
+
+ na = length(s_a);
+ nb = length(s_b);
+
+ if (na != nb)
+ error("Length of a and b differ (%i != %i)", na, nb);
+
+ ALLOC_REAL_VECTOR(na, s_ret, ret);
+
+ for (i = 0; i < na; ++i) {
+ double ca = dnorm(a[i], mean, sd, FALSE);
+ double cb = dnorm(b[i], mean, sd, FALSE);
+ double Ca = pnorm(a[i], mean, sd, FALSE, FALSE);
+ double Cb = pnorm(b[i], mean, sd, FALSE, FALSE);
+
+ ret[i] = mean + sd * ((ca - cb) / (Cb - Ca));
+ }
+ UNPROTECT(5);
+ return s_ret;
+}
More information about the Desire-commits
mailing list