[Desire-commits] r12 - in packages/truncnorm: . R man src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon May 19 10:38:15 CEST 2008
Author: olafm
Date: 2008-05-19 10:38:15 +0200 (Mon, 19 May 2008)
New Revision: 12
Modified:
packages/truncnorm/DESCRIPTION
packages/truncnorm/R/truncnorm.R
packages/truncnorm/man/dtruncnorm.Rd
packages/truncnorm/src/truncnorm.c
Log:
* Add vtruncnomr
Modified: packages/truncnorm/DESCRIPTION
===================================================================
--- packages/truncnorm/DESCRIPTION 2008-05-19 08:18:52 UTC (rev 11)
+++ packages/truncnorm/DESCRIPTION 2008-05-19 08:38:15 UTC (rev 12)
@@ -1,12 +1,12 @@
Package: truncnorm
-Version: 0.9.0
-Date: 2008-05-12
+Version: 0.9.1
+Date: 2008-05-19
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>
+Maintainer: Olaf Mersmann <olafm at statistik.uni-dortmund.de>
Depends: R (>= 2.7.0)
Description:
License: GPLv2
Modified: packages/truncnorm/R/truncnorm.R
===================================================================
--- packages/truncnorm/R/truncnorm.R 2008-05-19 08:18:52 UTC (rev 11)
+++ packages/truncnorm/R/truncnorm.R 2008-05-19 08:38:15 UTC (rev 12)
@@ -40,7 +40,17 @@
rtruncnorm <- function(n, a, b, mean=0, sd=1)
.External("rtruncnorm", n, a, b, mean, sd)
-etruncnorm <- function(a, b, mean=0, sd=1)
+etruncnorm <- function(a, b, mean, sd)
.External("etruncnorm", a, b, mean, sd)
+vtruncnorm <- function(a, b, mean, sd) {
+ na <- length(a)
+ if (na > 1) {
+ if (length(mean) == 1)
+ mean <- rep(mean, na)
+ if (length(sd) == 1)
+ sd <- rep(sd, na)
+ }
+ .External("etruncnorm", a, b, mean, sd)
+}
Modified: packages/truncnorm/man/dtruncnorm.Rd
===================================================================
--- packages/truncnorm/man/dtruncnorm.Rd 2008-05-19 08:18:52 UTC (rev 11)
+++ packages/truncnorm/man/dtruncnorm.Rd 2008-05-19 08:38:15 UTC (rev 12)
@@ -15,7 +15,8 @@
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)
+etruncnorm(a, b, mean, sd)
+etruncnorm(a, b, mean, sd)
}
\arguments{
\item{x,q}{vector of quantiles.}
@@ -33,8 +34,8 @@
\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.
+ generates random deviates, 'etruncnorm' gives the expected value and
+ 'vtruncnorm' the variance of the distirbution.
}
\references{ FIXME: Erstes auftauchen? }
\author{
Modified: packages/truncnorm/src/truncnorm.c
===================================================================
--- packages/truncnorm/src/truncnorm.c 2008-05-19 08:18:52 UTC (rev 11)
+++ packages/truncnorm/src/truncnorm.c 2008-05-19 08:38:15 UTC (rev 12)
@@ -51,7 +51,7 @@
c2 = sd * pnorm(b, mean, sd, FALSE, FALSE) - pnorm(a, mean, sd, FALSE, FALSE);
for (i = 0; i < n; ++i) {
- double t = x[i];
+ const double t = x[i];
if (a <= t && t <= b) { /* In range: */
c1 = dnorm(t, mean, sd, FALSE);
ret[i] = - c1 / c2;
@@ -67,7 +67,7 @@
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;
+ double *q, a, b, mean, sd, ca, cb, *ret;
UNPACK_REAL_VECTOR(args, s_q, q);
UNPACK_REAL_VALUE(args, s_a, a);
@@ -87,7 +87,7 @@
} else if (q[i] > b) {
ret[i] = 1.0;
} else {
- cx = pnorm(q[i], mean, sd, FALSE, FALSE);
+ const double cx = pnorm(q[i], mean, sd, FALSE, FALSE);
ret[i] = (cx - ca) / (cb - ca);
}
}
@@ -125,30 +125,69 @@
SEXP etruncnorm(SEXP args) {
- R_len_t i, na, nb;
+ R_len_t i, n_a, n_b, n_mean, n_sd;
SEXP s_a, s_b, s_mean, s_sd, s_ret;
- double *a, *b, mean, sd, *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);
+ UNPACK_REAL_VECTOR(args, s_mean, mean);
+ UNPACK_REAL_VECTOR(args, s_sd, sd);
- na = length(s_a);
- nb = length(s_b);
+ n_a = length(s_a);
+ n_b = length(s_b);
+ n_mean = length(s_mean);
+ n_sd = length(s_sd);
+ if (n_a != n_b || n_b != n_mean || n_mean != n_sd)
+ error("Length of a, b, mean or sd differ.");
+
+ ALLOC_REAL_VECTOR(n_a, s_ret, ret);
- if (na != nb)
- error("Length of a and b differ (%i != %i)", na, nb);
+ for (i = 0; i < n_a; ++i) {
+ const double ca = dnorm(a[i], mean[i], sd[i], FALSE);
+ const double cb = dnorm(b[i], mean[i], sd[i], FALSE);
+ const double Ca = pnorm(a[i], mean[i], sd[i], FALSE, FALSE);
+ const double Cb = pnorm(b[i], mean[i], sd[i], FALSE, FALSE);
+
+ ret[i] = mean[i] + sd[i] * ((ca - cb) / (Cb - Ca));
+ }
+ UNPROTECT(5);
+ return s_ret;
+}
+
- ALLOC_REAL_VECTOR(na, s_ret, ret);
+SEXP vtruncnorm(SEXP args) {
+ R_len_t i, n_a, n_b, n_mean, n_sd;
+ SEXP s_a, s_b, s_mean, s_sd, s_ret;
+ double *a, *b, *mean, *sd, *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);
+ UNPACK_REAL_VECTOR(args, s_a, a);
+ UNPACK_REAL_VECTOR(args, s_b, b);
+ UNPACK_REAL_VECTOR(args, s_mean, mean);
+ UNPACK_REAL_VECTOR(args, s_sd, sd);
+
+ n_a = length(s_a);
+ n_b = length(s_b);
+ n_mean = length(s_mean);
+ n_sd = length(s_sd);
+ if (n_a != n_b || n_b != n_mean || n_mean != n_sd)
+ error("Length of a, b, mean or sd differ.");
- ret[i] = mean + sd * ((ca - cb) / (Cb - Ca));
+ ALLOC_REAL_VECTOR(n_a, s_ret, ret);
+
+ for (i = 0; i < n_a; ++i) {
+ const double am = (a[i] - mean[i])/sd[i];
+ const double bm = (b[i] - mean[i])/sd[i];
+ const double ca = dnorm(am, 0.0, 1.0, FALSE);
+ const double cb = dnorm(bm, 0.0, 1.0, FALSE);
+ const double Ca = pnorm(am, 0.0, 1.0, FALSE, FALSE);
+ const double Cb = pnorm(bm, 0.0, 1.0, FALSE, FALSE);
+ const double v = sd[i] * sd[i];
+
+ const double d = Cb - Ca;
+ const double m1 = (ca - cb)/d;
+ const double m2 = (am*ca - bm*cb)/d;
+ ret[i] = (1.0 + m2 + m1*m1)*v;
}
UNPROTECT(5);
return s_ret;
More information about the Desire-commits
mailing list