[med-svn] [r-cran-truncnorm] 07/13: New upstream version 1.0-7
Andreas Tille
tille at debian.org
Thu Nov 30 14:39:04 UTC 2017
This is an automated email from the git hooks/post-receive script.
tille pushed a commit to branch master
in repository r-cran-truncnorm.
commit 7c9d03593df9a6429bfacbfe0b4da97fa8730556
Author: Andreas Tille <tille at debian.org>
Date: Thu Nov 30 15:19:49 2017 +0100
New upstream version 1.0-7
---
DESCRIPTION | 18 +++
MD5 | 10 ++
NAMESPACE | 13 ++
R/truncnorm.R | 33 +++++
debian/README.test | 8 --
debian/changelog | 17 ---
debian/compat | 1 -
debian/control | 23 ----
debian/copyright | 29 ----
debian/docs | 3 -
debian/rules | 4 -
debian/source/format | 1 -
debian/tests/control | 3 -
debian/tests/run-unit-test | 15 ---
debian/watch | 2 -
man/dtruncnorm.Rd | 59 +++++++++
src/rtruncnorm.c | 205 +++++++++++++++++++++++++++++
src/sexp_macros.h | 69 ++++++++++
src/truncnorm.c | 322 +++++++++++++++++++++++++++++++++++++++++++++
src/zeroin.c | 186 ++++++++++++++++++++++++++
src/zeroin.h | 12 ++
tests/sanity_checks.R | 161 +++++++++++++++++++++++
22 files changed, 1088 insertions(+), 106 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
new file mode 100644
index 0000000..0a6ccd9
--- /dev/null
+++ b/DESCRIPTION
@@ -0,0 +1,18 @@
+Package: truncnorm
+Version: 1.0-7
+Title: Truncated normal distribution
+Description: r/d/p/q functions for the truncated 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>
+ and Björn Bornkamp <bornkamp at statistik.tu-dortmund.de>
+Maintainer: Olaf Mersmann <olafm at statistik.uni-dortmund.de>
+Depends: R (>= 2.15.0)
+License: GPL-2
+Encoding: UTF-8
+Date:
+Packaged: 2014-01-21 08:56:20 UTC; olafm
+NeedsCompilation: yes
+Repository: CRAN
+Date/Publication: 2014-01-21 10:17:26
diff --git a/MD5 b/MD5
new file mode 100644
index 0000000..1f5e19d
--- /dev/null
+++ b/MD5
@@ -0,0 +1,10 @@
+dccda7b1843ccd076fef5a1efac46bf4 *DESCRIPTION
+6a17538f7760bcc8d34a441f44a32380 *NAMESPACE
+b1ae7a18e2a37ef976476ad504b8b4d6 *R/truncnorm.R
+43aaf800b50d70e17f63510a0d9edaf4 *man/dtruncnorm.Rd
+8d31aac61eb579e87150d3766697b8ef *src/rtruncnorm.c
+7c8d34cbe1c7a0e2783663e7ffb9327b *src/sexp_macros.h
+6b9f0732d4d643190af651ffbc369bfc *src/truncnorm.c
+06773d795a2fed9b539f1c248d12437a *src/zeroin.c
+dd074a0f18f05dda52777a57fd481e51 *src/zeroin.h
+768fe3c44063a1f13340084441bf40ba *tests/sanity_checks.R
diff --git a/NAMESPACE b/NAMESPACE
new file mode 100644
index 0000000..b6a9e43
--- /dev/null
+++ b/NAMESPACE
@@ -0,0 +1,13 @@
+export(dtruncnorm)
+export(ptruncnorm)
+export(qtruncnorm)
+export(rtruncnorm)
+export(etruncnorm)
+export(vtruncnorm)
+
+useDynLib(truncnorm,do_dtruncnorm)
+useDynLib(truncnorm,do_ptruncnorm)
+useDynLib(truncnorm,do_qtruncnorm)
+useDynLib(truncnorm,do_rtruncnorm)
+useDynLib(truncnorm,do_etruncnorm)
+useDynLib(truncnorm,do_vtruncnorm)
diff --git a/R/truncnorm.R b/R/truncnorm.R
new file mode 100644
index 0000000..6b5a913
--- /dev/null
+++ b/R/truncnorm.R
@@ -0,0 +1,33 @@
+##
+## 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=-Inf, b=Inf, mean=0, sd=1)
+ .Call("do_dtruncnorm", x, a, b, mean, sd)
+
+ptruncnorm <- function(q, a=-Inf, b=Inf, mean=0, sd=1)
+ .Call("do_ptruncnorm", q, a, b, mean, sd)
+
+qtruncnorm <- function(p, a=-Inf, b=Inf, mean=0, sd=1)
+ .Call("do_qtruncnorm", p, a, b, mean, sd)
+
+rtruncnorm <- function(n, a=-Inf, b=Inf, mean=0, sd=1) {
+ if (length(n) > 1)
+ n <- length(n)
+ if (length(n) > 1)
+ n <- length(n)
+ else if (!is.numeric(n))
+ stop("non-numeric argument n.")
+ .Call("do_rtruncnorm", as.integer(n), a, b, mean, sd)
+}
+
+etruncnorm <- function(a=-Inf, b=Inf, mean=0, sd=1)
+ .Call("do_etruncnorm", a, b, mean, sd)
+
+vtruncnorm <- function(a=-Inf, b=Inf, mean=0, sd=1)
+ .Call("do_vtruncnorm", a, b, mean, sd)
diff --git a/debian/README.test b/debian/README.test
deleted file mode 100644
index 55a9142..0000000
--- a/debian/README.test
+++ /dev/null
@@ -1,8 +0,0 @@
-Notes on how this package can be tested.
-────────────────────────────────────────
-
-To run the unit tests provided by the package you can do
-
- sh run-unit-test
-
-in this directory.
diff --git a/debian/changelog b/debian/changelog
deleted file mode 100644
index 9009e45..0000000
--- a/debian/changelog
+++ /dev/null
@@ -1,17 +0,0 @@
-r-cran-truncnorm (1.0-7-2) unstable; urgency=medium
-
- * Fix license
- Closes: #815687
- * Canonical homepage for CRAN
- * cme fix dpkg-control
- * Convert to dh-r
- * debhelper 10
- * d/watch: version=4
-
- -- Andreas Tille <tille at debian.org> Mon, 12 Dec 2016 11:16:51 +0100
-
-r-cran-truncnorm (1.0-7-1) unstable; urgency=low
-
- * Initial release (closes: #752164)
-
- -- Andreas Tille <tille at debian.org> Fri, 20 Jun 2014 11:47:43 +0200
diff --git a/debian/compat b/debian/compat
deleted file mode 100644
index f599e28..0000000
--- a/debian/compat
+++ /dev/null
@@ -1 +0,0 @@
-10
diff --git a/debian/control b/debian/control
deleted file mode 100644
index 223ee53..0000000
--- a/debian/control
+++ /dev/null
@@ -1,23 +0,0 @@
-Source: r-cran-truncnorm
-Maintainer: Debian Med Packaging Team <debian-med-packaging at lists.alioth.debian.org>
-Uploaders: Andreas Tille <tille at debian.org>
-Section: gnu-r
-Priority: optional
-Build-Depends: debhelper (>= 10),
- dh-r,
- r-base-dev
-Standards-Version: 3.9.8
-Vcs-Browser: https://anonscm.debian.org/viewvc/debian-med/trunk/packages/R/r-cran-truncnorm/trunk/
-Vcs-Svn: svn://anonscm.debian.org/debian-med/trunk/packages/R/r-cran-truncnorm/trunk/
-Homepage: https://cran.r-project.org/package=truncnorm
-
-Package: r-cran-truncnorm
-Architecture: any
-Depends: ${shlibs:Depends},
- ${misc:Depends},
- ${R:Depends}
-Recommends: ${R:Recommends}
-Suggests: ${R:Suggests}
-Description: GNU R truncated normal distribution
- This GNU R package provides r/d/p/q functions for the truncated normal
- distribution.
diff --git a/debian/copyright b/debian/copyright
deleted file mode 100644
index dc58101..0000000
--- a/debian/copyright
+++ /dev/null
@@ -1,29 +0,0 @@
-Format: https://www.debian.org/doc/packaging-manuals/copyright-format/1.0/
-Upstream-Name: truncnorm
-Upstream-Contact: Olaf Mersmann <olafm at statistik.uni-dortmund.de>
-Source: https://cran.r-project.org/package=truncnorm
-
-Files: *
-Copyright: 2012-2014 Heike Trautmann, Detlef Steuer, Olaf Mersmann, Björn Bornkamp
-License: GPL-2
-
-Files: debian/*
-Copyright: 2014-2016 Andreas Tille <tille at debian.org>
-License: GPL-2
-
-License: GPL-2
- This program is free software: you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation, either version 2 of the License.
- .
- This program is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU General Public License for more details.
- .
- You should have received a copy of the GNU General Public License
- along with this program. If not, see <http://www.gnu.org/licenses/>.
- .
- On Debian systems, the complete text of the GNU Public
- License can be found in `/usr/share/common-licenses/GPL-2'.
-
diff --git a/debian/docs b/debian/docs
deleted file mode 100644
index 3adf0d6..0000000
--- a/debian/docs
+++ /dev/null
@@ -1,3 +0,0 @@
-debian/README.test
-debian/tests/run-unit-test
-tests
diff --git a/debian/rules b/debian/rules
deleted file mode 100755
index 68d9a36..0000000
--- a/debian/rules
+++ /dev/null
@@ -1,4 +0,0 @@
-#!/usr/bin/make -f
-
-%:
- dh $@ --buildsystem R
diff --git a/debian/source/format b/debian/source/format
deleted file mode 100644
index 163aaf8..0000000
--- a/debian/source/format
+++ /dev/null
@@ -1 +0,0 @@
-3.0 (quilt)
diff --git a/debian/tests/control b/debian/tests/control
deleted file mode 100644
index 25377fc..0000000
--- a/debian/tests/control
+++ /dev/null
@@ -1,3 +0,0 @@
-Tests: run-unit-test
-Depends: @, r-cran-runit
-Restrictions: allow-stderr
diff --git a/debian/tests/run-unit-test b/debian/tests/run-unit-test
deleted file mode 100644
index 92689b6..0000000
--- a/debian/tests/run-unit-test
+++ /dev/null
@@ -1,15 +0,0 @@
-#!/bin/sh -e
-
-oname=truncnorm
-pkg=r-cran-`echo $oname | tr '[A-Z]' '[a-z]'`
-
-if [ "$ADTTMP" = "" ] ; then
- ADTTMP=`mktemp -d /tmp/${pkg}-test.XXXXXX`
- trap "rm -rf $ADTTMP" 0 INT QUIT ABRT PIPE TERM
-fi
-cd $ADTTMP
-cp /usr/share/doc/${pkg}/tests/* $ADTTMP
-find . -name "*.gz" -exec gunzip \{\} \;
-for runtest in `ls *.R` ; do
- R --no-save < $runtest
-done
diff --git a/debian/watch b/debian/watch
deleted file mode 100644
index 8c6226b..0000000
--- a/debian/watch
+++ /dev/null
@@ -1,2 +0,0 @@
-version=4
-http://cran.r-project.org/src/contrib/truncnorm_([-\d.]*)\.tar\.gz
diff --git a/man/dtruncnorm.Rd b/man/dtruncnorm.Rd
new file mode 100644
index 0000000..31f9ad1
--- /dev/null
+++ b/man/dtruncnorm.Rd
@@ -0,0 +1,59 @@
+\name{truncnorm}
+\alias{dtruncnorm}
+\alias{ptruncnorm}
+\alias{qtruncnorm}
+\alias{rtruncnorm}
+\alias{etruncnorm}
+\alias{vtruncnorm}
+\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=-Inf, b=Inf, mean = 0, sd = 1)
+ptruncnorm(q, a=-Inf, b=Inf, mean = 0, sd = 1)
+qtruncnorm(p, a=-Inf, b=Inf, mean = 0, sd = 1)
+rtruncnorm(n, a=-Inf, b=Inf, mean = 0, sd = 1)
+etruncnorm(a=-Inf, b=Inf, mean=0, sd=1)
+vtruncnorm(a=-Inf, b=Inf, mean=0, sd=1)
+}
+\arguments{
+ \item{x,q}{vector of quantiles.}
+ \item{p}{vector of probabilites.}
+ \item{n}{number of observations. If 'length(n) > 1', the length is
+ taken to be the number required.}
+ \item{a}{vector of lower bounds. These may be \code{-Inf}}
+ \item{b}{vector of upper bounds. These may be \code{Inf}}
+ \item{mean}{vector of means.}
+ \item{sd}{vector of standard deviations.}
+}
+\details{
+ If \code{mean} or \code{sd} are not specified they assume the default values of
+ \code{0} and \code{1}, respectively. The values of \code{a}, \code{b},
+ \code{mean} and \code{sd} are recycled as needed.
+}
+\value{
+ 'dtruncnorm' gives the density, 'ptruncnorm' gives the distribution
+ function, 'qtruncnorm' gives the quantile function, 'rtruncnorm'
+ generates random deviates, 'etruncnorm' gives the expected value and
+ 'vtruncnorm' the variance of the distribution.
+}
+\references{
+ The accept-reject sampler follows the description given in
+
+ Geweke, J. (1991). \emph{Efficient simulation from the multivariate normal
+ and student-t distributions subject to linear constraints}. In
+ Computing Science and Statistics: Proceedings of the 23rd Symposium on
+ the Interface, Ed. E. Keramidas and S. Kaufman, pp. 571-8. Fairfax
+ Station, VA: Interface Foundation of North America.
+}
+\author{
+ Heike Trautmann \email{trautmann at statistik.tu-dortmund.de},
+ Detlef Steuer \email{steuer at hsu-hamburg.de},
+ Olaf Mersmann \email{olafm at statistik.tu-dortmund.de} and
+ Björn Bornkamp \email{bornkamp at statistik.tu-dortmund.de} who donated a
+ much improved \code{rtruncnorm} implementation using an accept-reject sampler.
+}
+\keyword{distribution}
diff --git a/src/rtruncnorm.c b/src/rtruncnorm.c
new file mode 100644
index 0000000..b682623
--- /dev/null
+++ b/src/rtruncnorm.c
@@ -0,0 +1,205 @@
+/*
+ * rtruncnorm.c - Random truncated normal number generator.
+ *
+ * Authors:
+ * Björn Bornkamp <bornkamp at statistik.tu-dortmund.de>
+ * Olaf Mersmann <olafm at statistik.uni-dortmund.de>
+ */
+
+#include <R.h>
+#include <Rdefines.h>
+#include <Rmath.h>
+#include <Rinternals.h>
+#include <R_ext/Applic.h>
+#include <float.h>
+
+#include "sexp_macros.h"
+
+#define ALLOC_REAL_VECTOR(S, D, N) \
+ SEXP S; \
+ PROTECT(S = allocVector(REALSXP, N)); \
+ double *D = REAL(S);
+
+#ifndef MAX
+#define MAX(A, B) ((A>B)?(A):(B))
+#endif
+
+#include <R.h>
+#include <Rdefines.h>
+#include <Rmath.h>
+#include <Rinternals.h>
+#include <R_ext/Lapack.h>
+#include <R_ext/BLAS.h>
+#include <R_ext/Applic.h>
+
+#ifdef DEBUG
+#define SAMPLER_DEBUG(N, A, B) Rprintf("%8s(%f, %f)\n", N, A, B)
+#else
+#define SAMPLER_DEBUG(N, A, B)
+#endif
+
+static const double t1 = 0.15;
+static const double t2 = 2.18;
+static const double t3 = 0.725;
+static const double t4 = 0.45;
+
+/* Exponential rejection sampling (a,inf) */
+static R_INLINE double ers_a_inf(double a) {
+ SAMPLER_DEBUG("ers_a_inf", a, R_PosInf);
+ const double ainv = 1.0 / a;
+ double x, rho;
+ do{
+ x = rexp(ainv) + a; /* rexp works with 1/lambda */
+ rho = exp(-0.5 * pow((x - a), 2));
+ } while (runif(0, 1) > rho);
+ return x;
+}
+
+/* Exponential rejection sampling (a,b) */
+static R_INLINE double ers_a_b(double a, double b) {
+ SAMPLER_DEBUG("ers_a_b", a, b);
+ const double ainv = 1.0 / a;
+ double x, rho;
+ do{
+ x = rexp(ainv) + a; /* rexp works with 1/lambda */
+ rho = exp(-0.5 * pow((x-a), 2));
+ } while (runif(0, 1) > rho || x > b);
+ return x;
+}
+
+/* Normal rejection sampling (a,b) */
+static R_INLINE double nrs_a_b(double a, double b){
+ SAMPLER_DEBUG("nrs_a_b", a, b);
+ double x = -DBL_MAX;
+ while(x < a || x > b){
+ x = rnorm(0, 1);
+ }
+ return x;
+}
+
+/* Normal rejection sampling (a,inf) */
+static R_INLINE double nrs_a_inf(double a){
+ SAMPLER_DEBUG("nrs_a_inf", a, R_PosInf);
+ double x = -DBL_MAX;
+ while(x < a){
+ x = rnorm(0, 1);
+ }
+ return x;
+}
+
+/* Half-normal rejection sampling */
+double hnrs_a_b(double a, double b){
+ SAMPLER_DEBUG("hnrs_a_b", a, b);
+ double x = a - 1.0;
+ while(x < a || x > b) {
+ x = rnorm(0, 1);
+ x = fabs(x);
+ }
+ return x;
+}
+
+/* Uniform rejection sampling */
+static R_INLINE double urs_a_b(double a, double b){
+ SAMPLER_DEBUG("urs_a_b", a, b);
+ const double phi_a = dnorm(a, 0.0, 1.0, FALSE);
+ double x = 0.0, u = 0.0;
+
+ /* Upper bound of normal density on [a, b] */
+ const double ub = a < 0 && b > 0 ? M_1_SQRT_2PI : phi_a;
+ do{
+ x = runif(a, b);
+ } while (runif(0, 1) * ub > dnorm(x,0,1,0));
+ return x;
+}
+
+/* Previously this was refered to as type 1 sampling: */
+static inline double r_lefttruncnorm(double a, double mean, double sd) {
+ const double alpha = (a - mean) / sd;
+ if (alpha < t4) {
+ return mean + sd * nrs_a_inf(alpha);
+ } else {
+ return mean + sd * ers_a_inf(alpha);
+ }
+}
+
+static R_INLINE double r_righttruncnorm(double b, double mean, double sd) {
+ const double beta = (b - mean) / sd;
+ /* Exploit symmetry: */
+ return mean - sd * r_lefttruncnorm(-beta, 0.0, 1.0);
+}
+
+static R_INLINE double r_truncnorm(double a, double b, double mean, double sd) {
+ const double alpha = (a - mean) / sd;
+ const double beta = (b - mean) / sd;
+ const double phi_a = dnorm(alpha, 0.0, 1.0, FALSE);
+ const double phi_b = dnorm(beta, 0.0, 1.0, FALSE);
+ if (beta <= alpha) {
+ return NA_REAL;
+ } else if (alpha <= 0 && 0 <= beta) { /* 2 */
+ if (phi_a <= t1 || phi_b <= t1) { /* 2 (a) */
+ return mean + sd * nrs_a_b(alpha, beta);
+ } else { /* 2 (b) */
+ return mean + sd * urs_a_b(alpha, beta);
+ }
+ } else if (alpha > 0) { /* 3 */
+ if (phi_a / phi_b <= t2) { /* 3 (a) */
+ return mean + sd * urs_a_b(alpha, beta);
+ } else {
+ if (alpha < t3) { /* 3 (b) */
+ return mean + sd * hnrs_a_b(alpha, beta);
+ } else { /* 3 (c) */
+ return mean + sd * ers_a_b(alpha, beta);
+ }
+ }
+ } else { /* 3s */
+ if (phi_b / phi_a <= t2) { /* 3s (a) */
+ return mean - sd * urs_a_b(-beta, -alpha);
+ } else {
+ if (beta > -t3) { /* 3s (b) */
+ return mean - sd * hnrs_a_b(-beta, -alpha);
+ } else { /* 3s (c) */
+ return mean - sd * ers_a_b(-beta, -alpha);
+ }
+ }
+ }
+}
+
+SEXP do_rtruncnorm(SEXP s_n, SEXP s_a, SEXP s_b, SEXP s_mean, SEXP s_sd) {
+ R_len_t i, nn;
+ UNPACK_INT(s_n, n);
+ if (NA_INTEGER == n)
+ error("n is NA - aborting.");
+ UNPACK_REAL_VECTOR(s_a , a , n_a);
+ UNPACK_REAL_VECTOR(s_b , b , n_b);
+ UNPACK_REAL_VECTOR(s_mean, mean, n_mean);
+ UNPACK_REAL_VECTOR(s_sd , sd , n_sd);
+
+ nn = MAX(n, MAX(MAX(n_a, n_b), MAX(n_mean, n_sd)));
+ ALLOC_REAL_VECTOR(s_ret, ret, nn);
+
+ GetRNGstate();
+ for (i = 0; i < nn; ++i) {
+ const double ca = a[i % n_a];
+ const double cb = b[i % n_b];
+ const double cmean = mean[i % n_mean];
+ const double csd = sd[i % n_sd];
+
+ if (R_FINITE(ca) && R_FINITE(cb)) {
+ ret[i] = r_truncnorm(ca, cb, cmean, csd);
+ } else if (R_NegInf == ca && R_FINITE(cb)) {
+ ret[i] = r_righttruncnorm(cb, cmean, csd);
+ } else if (R_FINITE(ca) && R_PosInf == cb) {
+ ret[i] = r_lefttruncnorm(ca, cmean, csd);
+ } else if (R_NegInf == ca && R_PosInf == cb) {
+ ret[i] = rnorm(cmean, csd);
+ } else {
+ ret[i] = NA_REAL;
+ }
+ R_CheckUserInterrupt();
+ }
+ PutRNGstate();
+ UNPROTECT(1); /* s_ret */
+ return s_ret;
+}
+
+
diff --git a/src/sexp_macros.h b/src/sexp_macros.h
new file mode 100644
index 0000000..3a12d1c
--- /dev/null
+++ b/src/sexp_macros.h
@@ -0,0 +1,69 @@
+/*
+ * sexp_macros.h - helper macros for SEXPs
+ *
+ * Collection of useful macros to handle S expressions. Most of these
+ * are used to unpack arguments passed in via the .Call() or
+ * .External() interface.
+ *
+ * Author:
+ * Olaf Mersmann (OME) <olafm at statistik.tu-dortmund.de>
+ */
+
+#if !defined(__SEXP_MACROS_H__)
+#define __SEXP_MACROS_H__
+
+#include <R.h>
+#include <Rinternals.h>
+
+#define CHECK_ARG_IS_REAL_MATRIX(A) \
+ if (!isReal(A) || !isMatrix(A)) \
+ error("Argument '" #A "' is not a real matrix.");
+
+#define CHECK_ARG_IS_REAL_VECTOR(A) \
+ if (!isReal(A) || !isVector(A)) \
+ error("Argument '" #A "' is not a real vector.");
+
+#define CHECK_ARG_IS_INT_VECTOR(A) \
+ if (!isInteger(A) || !isVector(A)) \
+ error("Argument '" #A "' is not an integer vector.");
+
+/*
+ * Unpack a real matrix stored in SEXP S.
+ */
+#define UNPACK_REAL_MATRIX(S, D, N, K) \
+ CHECK_ARG_IS_REAL_MATRIX(S); \
+ double *D = REAL(S); \
+ const R_len_t N = nrows(S); \
+ const R_len_t K = ncols(S);
+
+/*
+ * Unpack a real vector stored in SEXP S.
+ */
+#define UNPACK_REAL_VECTOR(S, D, N) \
+ CHECK_ARG_IS_REAL_VECTOR(S); \
+ double *D = REAL(S); \
+ const R_len_t N = length(S);
+
+/*
+ * Unpack a single real stored in SEXP S.
+ */
+#define UNPACK_REAL(S, D) \
+ CHECK_ARG_IS_REAL_VECTOR(S); \
+ double D = REAL(S)[0]; \
+
+/*
+ * Unpack an integer vector stored in SEXP S.
+ */
+#define UNPACK_INT_VECTOR(S, I, N) \
+ CHECK_ARG_IS_INT_VECTOR(S); \
+ int *I = INTEGER(S); \
+ const R_len_t N = length(S);
+
+/*
+ * Unpack a single integer stored in SEXP S.
+ */
+#define UNPACK_INT(S, I) \
+ CHECK_ARG_IS_INT_VECTOR(S); \
+ int I = INTEGER(S)[0]; \
+
+#endif
diff --git a/src/truncnorm.c b/src/truncnorm.c
new file mode 100644
index 0000000..52650f7
--- /dev/null
+++ b/src/truncnorm.c
@@ -0,0 +1,322 @@
+/*
+ * 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>
+
+#include "sexp_macros.h"
+#include "zeroin.h"
+
+#define ALLOC_REAL_VECTOR(S, D, N) \
+ SEXP S; \
+ PROTECT(S = allocVector(REALSXP, N)); \
+ double *D = REAL(S);
+
+#ifndef MAX
+#define MAX(A, B) ((A>B)?(A):(B))
+#endif
+
+/*
+ * These routines calculate the expected value and variance of the
+ * left, right and doubly truncated normal distribution. The only
+ * tricky bit is the calculation of the variance of the doubly
+ * truncated normal distribution. We use a decompostion of the
+ * variance of a mixture of distributions to here for numerical
+ * reasons. For details see:
+ *
+ * Foulley JL. A completion simulator for the two-sided truncated
+ * normal distribution. Genetics, selection, evolution 2000
+ * Nov-Dec;32(6): p. 631-635.
+ */
+static R_INLINE double e_lefttruncnorm(double a, double mean, double sd) {
+ const double alpha = (a - mean) / sd;
+ const double phi_a = dnorm(alpha, 0.0, 1.0, FALSE);
+ const double Phi_a = pnorm(alpha, 0.0, 1.0, TRUE, FALSE);
+ double res = mean + sd * (phi_a / (1.0 - Phi_a));
+ return res;
+}
+
+static R_INLINE double e_truncnorm(double a, double b, double mean, double sd) {
+ double delta_phi = 0.0, delta_Phi = 0.0;
+
+ const double alpha = (a - mean) / sd;
+ const double beta = (b - mean) / sd;
+ const double phi_a = dnorm(alpha, 0.0, 1.0, TRUE);
+ const double Phi_a = pnorm(alpha, 0.0, 1.0, TRUE, TRUE);
+ const double phi_b = dnorm(beta, 0.0, 1.0, TRUE);
+ const double Phi_b = pnorm(beta, 0.0, 1.0, TRUE, TRUE);
+
+ if (phi_b < phi_a) {
+ delta_phi = logspace_sub(phi_a, phi_b);
+ } else {
+ sd = -sd;
+ delta_phi = logspace_sub(phi_b, phi_a);
+ }
+
+ if (Phi_b > Phi_a) {
+ sd = -sd;
+ delta_Phi = logspace_sub(Phi_b, Phi_a);
+ } else {
+ delta_Phi = logspace_sub(Phi_a, Phi_b);
+ }
+ /* Rprintf("pb - pa = dp: %.16f - %.16f = %.16f\n", phi_b, phi_a, delta_phi); */
+ /* Rprintf("Pa - Pb = dP: %.16f - %.16f = %.16f\n", Phi_a, Phi_b, delta_Phi); */
+
+ return mean + sd *-exp(delta_phi - delta_Phi);
+}
+
+static R_INLINE double e_righttruncnorm(double b, double mean, double sd) {
+ const double beta = (b - mean) / sd;
+ const double phi_b = dnorm(beta, 0.0, 1.0, FALSE);
+ const double Phi_b = pnorm(beta, 0.0, 1.0, TRUE, FALSE);
+ return mean + sd * (-phi_b / Phi_b);
+}
+
+static R_INLINE double v_lefttruncnorm(double a, double mean, double sd) {
+ const double alpha = (a - mean) / sd;
+ const double phi_a = dnorm(alpha, 0.0, 1.0, FALSE);
+ const double Phi_a = pnorm(alpha, 0.0, 1.0, TRUE, FALSE);
+ const double lambda = phi_a / (1.0 - Phi_a);
+ return (sd*sd*(1.0 - lambda * (lambda - alpha)));
+}
+
+static R_INLINE double v_righttruncnorm(double b, double mean, double sd) {
+ return (v_lefttruncnorm(-b, -mean, sd));
+}
+
+static R_INLINE double v_truncnorm(double a, double b, double mean, double sd) {
+ const double pi1 = pnorm(a, mean, sd, TRUE, FALSE);
+ const double pi2 = pnorm(b, mean, sd, TRUE, FALSE) - pnorm(a, mean, sd, TRUE, FALSE);
+ const double pi3 = pnorm(b, mean, sd, FALSE, FALSE); /* 1 - F(b) */
+
+ const double e1 = e_righttruncnorm(a, mean, sd);
+ const double e2 = e_truncnorm(a, b, mean, sd);
+ const double e3 = e_lefttruncnorm(b, mean, sd);
+
+ const double v = sd * sd;
+ const double v1 = v_righttruncnorm(a, mean, sd);
+ const double v3 = v_lefttruncnorm(b, mean, sd);
+
+ const double c1 = pi1 * (v1 + (e1 - mean)*(e1 - mean));
+ const double c3 = pi3 * (v3 + (e3 - mean)*(e3 - mean));
+
+ return (v - c1 - c3) / pi2 - (e2 - mean)*(e2 - mean);
+}
+
+static R_INLINE double ptruncnorm(const double q,
+ const double a, const double b,
+ const double mean, const double sd) {
+ if (q < a) {
+ return 0.0;
+ } else if (q > b) {
+ return 1.0;
+ } else {
+ const double c1 = pnorm(q, mean, sd, TRUE, FALSE);
+ const double c2 = pnorm(a, mean, sd, TRUE, FALSE);
+ const double c3 = pnorm(b, mean, sd, TRUE, FALSE);
+ return (c1 - c2) / (c3 - c2);
+ }
+}
+
+typedef struct {
+ double a, b, mean, sd, p;
+} qtn;
+
+/* qtmin - helper function to calculate quantiles of the truncated
+ * normal distribution.
+ *
+ * The root of this function is the desired quantile, given that *p
+ * defines a truncated normal distribution and the desired
+ * quantile. This function increases monotonically in x and is
+ * positive for x=a and negative for x=b if 0 < p < 1.
+ */
+double qtmin(double x, void *p) {
+ qtn *t = (qtn *)p;
+ return ptruncnorm(x, t->a, t->b, t->mean, t->sd) - t->p;
+}
+
+SEXP do_dtruncnorm(SEXP s_x, SEXP s_a, SEXP s_b, SEXP s_mean, SEXP s_sd) {
+ R_len_t i, n;
+ UNPACK_REAL_VECTOR(s_x , x , n_x);
+ UNPACK_REAL_VECTOR(s_a , a , n_a);
+ UNPACK_REAL_VECTOR(s_b , b , n_b);
+ UNPACK_REAL_VECTOR(s_mean, mean, n_mean);
+ UNPACK_REAL_VECTOR(s_sd , sd , n_sd);
+
+ n = MAX(MAX(MAX(n_x, n_a), MAX(n_b, n_mean)), n_sd);
+ ALLOC_REAL_VECTOR(s_ret, ret, n);
+
+ for (i = 0; i < n; ++i) {
+ const double ca = a[i % n_a];
+ const double cb = b[i % n_b];
+ const double cx = x[i];
+ if (ca <= cx && cx <= cb) { /* In range: */
+ const double cmean = mean[i % n_mean];
+ const double csd = sd[i % n_sd];
+
+ const double c1 = pnorm(ca, cmean, csd, TRUE, FALSE);
+ const double c2 = pnorm(cb, cmean, csd, TRUE, FALSE);
+ const double c3 = csd * (c2 - c1);
+ const double c4 = dnorm((cx-cmean)/csd, 0.0, 1.0, FALSE);
+ ret[i] = c4 / c3;
+ } else { /* Truncated: */
+ ret[i] = 0.0;
+ }
+ R_CheckUserInterrupt();
+ }
+ UNPROTECT(1); /* s_ret */
+ return s_ret;
+}
+
+SEXP do_ptruncnorm(SEXP s_q, SEXP s_a, SEXP s_b, SEXP s_mean, SEXP s_sd) {
+ R_len_t i, n;
+ UNPACK_REAL_VECTOR(s_q , q , n_q);
+ UNPACK_REAL_VECTOR(s_a , a , n_a);
+ UNPACK_REAL_VECTOR(s_b , b , n_b);
+ UNPACK_REAL_VECTOR(s_mean, mean, n_mean);
+ UNPACK_REAL_VECTOR(s_sd , sd , n_sd);
+
+ n = MAX(MAX(MAX(n_q, n_a), MAX(n_b, n_mean)), n_sd);
+ ALLOC_REAL_VECTOR(s_ret, ret, n);
+
+ for (i = 0; i < n; ++i) {
+ const double cq = q[i % n_q];
+ const double ca = a[i % n_a];
+ const double cb = b[i % n_b];
+ const double cmean = mean[i % n_mean];
+ const double csd = sd[i % n_sd];
+ ret[i] = ptruncnorm(cq, ca, cb, cmean, csd);
+ R_CheckUserInterrupt();
+ }
+ UNPROTECT(1); /* s_ret */
+ return s_ret;
+}
+
+SEXP do_qtruncnorm(SEXP s_p, SEXP s_a, SEXP s_b, SEXP s_mean, SEXP s_sd) {
+ R_len_t i, n;
+ qtn t;
+ double tol;
+ int maxit;
+ UNPACK_REAL_VECTOR(s_p , p , n_p);
+ UNPACK_REAL_VECTOR(s_a , a , n_a);
+ UNPACK_REAL_VECTOR(s_b , b , n_b);
+ UNPACK_REAL_VECTOR(s_mean, mean, n_mean);
+ UNPACK_REAL_VECTOR(s_sd , sd , n_sd);
+
+ n = MAX(MAX(MAX(n_p, n_a), MAX(n_b, n_mean)), n_sd);
+ ALLOC_REAL_VECTOR(s_ret, ret, n);
+
+ for (i = 0; i < n; ++i) {
+ const double cp = p[i % n_p];
+ const double ca = a[i % n_a];
+ const double cb = b[i % n_b];
+ const double cmean = mean[i % n_mean];
+ const double csd = sd[i % n_sd];
+ if (cp == 0.0) {
+ ret[i] = ca;
+ } else if (cp == 1.0) {
+ ret[i] = cb;
+ } else if (cp < 0.0 || cp > 1.0) {
+ ret[i] = R_NaN;
+ } else if (ca == R_NegInf && cb == R_PosInf) {
+ ret[i] = qnorm(cp, cmean, csd, TRUE, FALSE);
+ } else {
+ /* We need to possible adjust ca and cb for R_zeroin(),
+ * because R_zeroin() requires finite bounds and ca or cb (but
+ * not both, see above) may be infinite. In that case, we use
+ * a simple stepping out procedure to find a lower or upper
+ * bound from which to begin the search.
+ */
+ double lower = ca, upper = cb;
+ if (lower == R_NegInf) {
+ lower = -1;
+ while(ptruncnorm(lower, ca, cb, cmean, csd) - cp >= 0) lower *= 2.0;
+ } else if (upper == R_PosInf) {
+ upper = 1;
+ while(ptruncnorm(upper, ca, cb, cmean, csd) - cp <= 0) upper *= 2.0;
+ }
+ t.a = ca; t.b = cb; t.mean = cmean; t.sd = csd; t.p = cp; maxit = 200;
+ tol = 0.0; /* Set tolerance! */
+ ret[i] = truncnorm_zeroin(lower, upper,
+ qtmin(lower, &t), qtmin(upper, &t),
+ &qtmin, &t, &tol, &maxit);
+ }
+ R_CheckUserInterrupt();
+ }
+ UNPROTECT(1); /* s_ret */
+ return s_ret;
+}
+
+SEXP do_etruncnorm(SEXP s_a, SEXP s_b, SEXP s_mean, SEXP s_sd) {
+ R_len_t i, n;
+ UNPACK_REAL_VECTOR(s_a , a , n_a);
+ UNPACK_REAL_VECTOR(s_b , b , n_b);
+ UNPACK_REAL_VECTOR(s_mean, mean, n_mean);
+ UNPACK_REAL_VECTOR(s_sd , sd , n_sd);
+
+ n = MAX(MAX(n_a, n_b), MAX(n_mean, n_sd));
+ ALLOC_REAL_VECTOR(s_ret, ret, n);
+
+ for (i = 0; i < n; ++i) {
+ const double ca = a[i % n_a];
+ const double cb = b[i % n_b];
+ const double cmean = mean[i % n_mean];
+ const double csd = sd[i % n_sd];
+
+ if (R_FINITE(ca) && R_FINITE(cb)) {
+ ret[i] = e_truncnorm(ca, cb, cmean, csd);
+ } else if (R_NegInf == ca && R_FINITE(cb)) {
+ ret[i] = e_righttruncnorm(cb, cmean, csd);
+ } else if (R_FINITE(ca) && R_PosInf == cb) {
+ ret[i] = e_lefttruncnorm(ca, cmean, csd);
+ } else if (R_NegInf == ca && R_PosInf == cb) {
+ ret[i] = cmean;
+ } else {
+ ret[i] = NA_REAL;
+ }
+ R_CheckUserInterrupt();
+ }
+ UNPROTECT(1); /* s_ret */
+ return s_ret;
+}
+
+SEXP do_vtruncnorm(SEXP s_a, SEXP s_b, SEXP s_mean, SEXP s_sd) {
+ R_len_t i, n;
+ UNPACK_REAL_VECTOR(s_a , a , n_a);
+ UNPACK_REAL_VECTOR(s_b , b , n_b);
+ UNPACK_REAL_VECTOR(s_mean, mean, n_mean);
+ UNPACK_REAL_VECTOR(s_sd , sd , n_sd);
+
+ n = MAX(MAX(n_a, n_b), MAX(n_mean, n_sd));
+ ALLOC_REAL_VECTOR(s_ret, ret, n);
+
+ for (i = 0; i < n; ++i) {
+ const double ca = a[i % n_a];
+ const double cb = b[i % n_b];
+ const double cmean = mean[i % n_mean];
+ const double csd = sd[i % n_sd];
+
+ if (R_FINITE(ca) && R_FINITE(cb)) {
+ ret[i] = v_truncnorm(ca, cb, cmean, csd);
+ } else if (R_NegInf == ca && R_FINITE(cb)) {
+ ret[i] = v_righttruncnorm(cb, cmean, csd);
+ } else if (R_FINITE(ca) && R_PosInf == cb) {
+ ret[i] = v_lefttruncnorm(ca, cmean, csd);
+ } else if (R_NegInf == ca && R_PosInf == cb) {
+ ret[i] = csd * csd;
+ } else {
+ ret[i] = NA_REAL;
+ }
+ R_CheckUserInterrupt();
+ }
+ UNPROTECT(1); /* s_ret */
+ return s_ret;
+}
diff --git a/src/zeroin.c b/src/zeroin.c
new file mode 100644
index 0000000..e4d3ded
--- /dev/null
+++ b/src/zeroin.c
@@ -0,0 +1,186 @@
+/* This code was taken from the main R distribution.
+ */
+
+/*
+ * R : A Computer Language for Statistical Data Analysis
+ * Copyright (C) 1999, 2001 the R Core Team
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, a copy is available at
+ * http://www.r-project.org/Licenses/
+ */
+
+/* from NETLIB c/brent.shar with max.iter, add'l info and convergence
+ details hacked in by Peter Dalgaard */
+
+/*************************************************************************
+ * C math library
+ * function ZEROIN - obtain a function zero within the given range
+ *
+ * Input
+ * double zeroin(ax,bx,f,info,Tol,Maxit)
+ * double ax; Root will be seeked for within
+ * double bx; a range [ax,bx]
+ * double (*f)(double x, void *info); Name of the function whose zero
+ * will be seeked for
+ * void *info; Add'l info passed to f
+ * double *Tol; Acceptable tolerance for the root
+ * value.
+ * May be specified as 0.0 to cause
+ * the program to find the root as
+ * accurate as possible
+ *
+ * int *Maxit; Max. iterations
+ *
+ *
+ * Output
+ * Zeroin returns an estimate for the root with accuracy
+ * 4*EPSILON*abs(x) + tol
+ * *Tol returns estimated precision
+ * *Maxit returns actual # of iterations, or -1 if maxit was
+ * reached without convergence.
+ *
+ * Algorithm
+ * G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical
+ * computations. M., Mir, 1980, p.180 of the Russian edition
+ *
+ * The function makes use of the bisection procedure combined with
+ * the linear or quadric inverse interpolation.
+ * At every step program operates on three abscissae - a, b, and c.
+ * b - the last and the best approximation to the root
+ * a - the last but one approximation
+ * c - the last but one or even earlier approximation than a that
+ * 1) |f(b)| <= |f(c)|
+ * 2) f(b) and f(c) have opposite signs, i.e. b and c confine
+ * the root
+ * At every step Zeroin selects one of the two new approximations, the
+ * former being obtained by the bisection procedure and the latter
+ * resulting in the interpolation (if a,b, and c are all different
+ * the quadric interpolation is utilized, otherwise the linear one).
+ * If the latter (i.e. obtained by the interpolation) point is
+ * reasonable (i.e. lies within the current interval [b,c] not being
+ * too close to the boundaries) it is accepted. The bisection result
+ * is used in the other case. Therefore, the range of uncertainty is
+ * ensured to be reduced at least by the factor 1.6
+ */
+
+#include <float.h>
+#include <math.h>
+
+#define EPSILON DBL_EPSILON
+
+double truncnorm_zeroin( /* An estimate of the root */
+ double ax, /* Left border | of the range */
+ double bx, /* Right border| the root is seeked*/
+ double fa, double fb, /* f(a), f(b) */
+ double (*f)(double x, void *info), /* Function under investigation */
+ void *info, /* Add'l info passed on to f */
+ double *Tol, /* Acceptable tolerance */
+ int *Maxit) /* Max # of iterations */
+{
+ double a, b, c, fc; /* Abscissae, descr. see above, f(c) */
+ double tol;
+ int maxit;
+
+ a = ax; b = bx;
+ c = a; fc = fa;
+ maxit = *Maxit + 1; tol = * Tol;
+
+ /* First test if we have found a root at an endpoint */
+ if(fa == 0.0) {
+ *Tol = 0.0;
+ *Maxit = 0;
+ return a;
+ }
+ if(fb == 0.0) {
+ *Tol = 0.0;
+ *Maxit = 0;
+ return b;
+ }
+
+ while(maxit--) /* Main iteration loop */
+ {
+ double prev_step = b-a; /* Distance from the last but one
+ to the last approximation */
+ double tol_act; /* Actual tolerance */
+ double p; /* Interpolation step is calcu- */
+ double q; /* lated in the form p/q; divi-
+ * sion operations is delayed
+ * until the last moment */
+ double new_step; /* Step at this iteration */
+
+ if( fabs(fc) < fabs(fb) )
+ { /* Swap data for b to be the */
+ a = b; b = c; c = a; /* best approximation */
+ fa=fb; fb=fc; fc=fa;
+ }
+ tol_act = 2*EPSILON*fabs(b) + tol/2;
+ new_step = (c-b)/2;
+
+ if( fabs(new_step) <= tol_act || fb == (double)0 )
+ {
+ *Maxit -= maxit;
+ *Tol = fabs(c-b);
+ return b; /* Acceptable approx. is found */
+ }
+
+ /* Decide if the interpolation can be tried */
+ if( fabs(prev_step) >= tol_act /* If prev_step was large enough*/
+ && fabs(fa) > fabs(fb) ) { /* and was in true direction,
+ * Interpolation may be tried */
+ register double t1,cb,t2;
+ cb = c-b;
+ if( a==c ) { /* If we have only two distinct */
+ /* points linear interpolation */
+ t1 = fb/fa; /* can only be applied */
+ p = cb*t1;
+ q = 1.0 - t1;
+ }
+ else { /* Quadric inverse interpolation*/
+
+ q = fa/fc; t1 = fb/fc; t2 = fb/fa;
+ p = t2 * ( cb*q*(q-t1) - (b-a)*(t1-1.0) );
+ q = (q-1.0) * (t1-1.0) * (t2-1.0);
+ }
+ if( p>(double)0 ) /* p was calculated with the */
+ q = -q; /* opposite sign; make p positive */
+ else /* and assign possible minus to */
+ p = -p; /* q */
+
+ if( p < (0.75*cb*q-fabs(tol_act*q)/2) /* If b+p/q falls in [b,c]*/
+ && p < fabs(prev_step*q/2) ) /* and isn't too large */
+ new_step = p/q; /* it is accepted
+ * If p/q is too large then the
+ * bisection procedure can
+ * reduce [b,c] range to more
+ * extent */
+ }
+
+ if( fabs(new_step) < tol_act) { /* Adjust the step to be not less*/
+ if( new_step > (double)0 ) /* than tolerance */
+ new_step = tol_act;
+ else
+ new_step = -tol_act;
+ }
+ a = b; fa = fb; /* Save the previous approx. */
+ b += new_step; fb = (*f)(b, info); /* Do step to a new approxim. */
+ if( (fb > 0 && fc > 0) || (fb < 0 && fc < 0) ) {
+ /* Adjust c for it to have a sign opposite to that of b */
+ c = a; fc = fa;
+ }
+ }
+ /* failed! */
+ *Tol = fabs(c-b);
+ *Maxit = -1;
+ return b;
+}
diff --git a/src/zeroin.h b/src/zeroin.h
new file mode 100644
index 0000000..4f9d65e
--- /dev/null
+++ b/src/zeroin.h
@@ -0,0 +1,12 @@
+#ifndef __zeroin_h__
+#define __zeroin_h__
+double truncnorm_zeroin( /* An estimate of the root */
+ double ax, /* Left border | of the range */
+ double bx, /* Right border| the root is seeked*/
+ double fa, double fb, /* f(a), f(b) */
+ double (*f)(double x, void *info), /* Function under investigation */
+ void *info, /* Add'l info passed on to f */
+ double *Tol, /* Acceptable tolerance */
+ int *Maxit);
+
+#endif
diff --git a/tests/sanity_checks.R b/tests/sanity_checks.R
new file mode 100644
index 0000000..cff9550
--- /dev/null
+++ b/tests/sanity_checks.R
@@ -0,0 +1,161 @@
+require("truncnorm")
+
+################################################################################
+## Check d/e/vtruncnorm all in one function:
+check_dev <- function(a, b, mean=0, sd=1) {
+ e <- etruncnorm(a, b, mean, sd)
+ v <- vtruncnorm(a, b, mean, sd)
+
+ ee <- integrate(function(x) x * dtruncnorm(x, a, b, mean, sd), a, b)$value
+ if (abs(e - ee)/e > 0.00005) {
+ message(sprintf("FAIL: etruncnorm(%4.1f, %4.1f, %4.1f, %4.1f) - mismatch %f vs. %f",
+ a, b, mean, sd, ee, e))
+ }
+ ev <- integrate(function(x) (x-ee)^2 * dtruncnorm(x, a, b, mean, sd), a, b)$value
+ if (abs(v - ev)/v > 0.00005) {
+ message(sprintf("FAIL: vtruncnorm(%4.1f, %4.1f, %4.1f, %4.1f) - mismatch %f vs. %f",
+ a, b, mean, sd, ev, v))
+ }
+}
+
+## Left truncated:
+check_dev(-3, Inf, 0, 1)
+check_dev(-2, Inf, 1, 1)
+check_dev( 2, Inf, 0, 1)
+check_dev( 3, Inf, 1, 1)
+
+check_dev(-3, Inf, 0, 2)
+check_dev(-2, Inf, 1, 2)
+check_dev( 2, Inf, 0, 2)
+check_dev( 3, Inf, 1, 2)
+
+## Doubly truncated:
+check_dev(-3.0, -2.5, 0, 1)
+check_dev(-3.0, -1.5, 0, 1)
+check_dev(-3.0, -0.5, 0, 1)
+check_dev(-3.0, 0.5, 0, 1)
+
+check_dev(0.0, 0.5, 0, 1)
+check_dev(0.0, 1.5, 0, 1)
+check_dev(0.0, 2.5, 0, 1)
+check_dev(0.0, 3.5, 0, 1)
+
+## Extreme cases:
+check_dev( 0.0, 1.0, 0.0, 10.0)
+check_dev( 0.0, 1.0, 5.0, 1.0)
+check_dev(-1.0, 0.0, 0.0, 10.0)
+check_dev( 0.0, 1.0, -5.0, 1.0)
+
+################################################################################
+## Sanity checks on random number generators
+check_r <- function(a, b, mean, sd, n=10000) {
+ x <- rtruncnorm(n, a, b, mean, sd)
+ if (!all(x > a & x < b)) {
+ message(sprintf("FAIL: rtruncnorm(%i, %4.1f, %4.1f, %4.1f, %4.1f) - bounds",
+ n, a, b, mean, sd))
+ }
+
+ ## Check to make sure mean and variance have the correct magnitude.
+ e.x <- mean(x)
+ e <- etruncnorm(a, b, mean, sd)
+ if (abs(e.x - e)/sd > 0.05) {
+ message(sprintf("FAIL: rtruncnorm(%i, %4.1f, %4.1f, %4.1f, %4.1f) - mean %f vs. %f",
+ n, a, b, mean, sd, e.x, e))
+ }
+ sd.x <- sd(x)
+ sd <- sqrt(vtruncnorm(a, b, mean, sd))
+ if (abs(sd.x - sd)/sd.x > 0.05) {
+ message(sprintf("FAIL: rtruncnorm(%i, %4.1f, %4.1f, %4.1f, %4.1f) - variance %f vs. %f",
+ n, a, b, mean, sd, sd.x, sd))
+ }
+}
+
+## rtruncnorm == rnorm:
+check_r(-Inf, Inf, 0, 1)
+
+## 0 in (a, b):
+check_r(-1, 1, 0, 1)
+check_r(-1, 1, 1, 1)
+check_r(-1, 1, 0, 2)
+
+## 0 < (a, b):
+check_r(1, 2, 0, 1)
+check_r(1, 2, 1, 1)
+check_r(1, 2, 0, 2)
+
+## 0 > (a, b):
+check_r(-2, -1, 0, 1)
+check_r(-2, -1, 1, 1)
+check_r(-2, -1, 0, 2)
+
+## left truncation:
+check_r(-2, Inf, 0, 1)
+check_r(-2, Inf, 1, 1)
+check_r(-2, Inf, 0, 2)
+check_r( 0, Inf, 0, 1)
+check_r( 0, Inf, 1, 1)
+check_r( 0, Inf, 0, 2)
+check_r( 2, Inf, 0, 1)
+check_r( 2, Inf, 1, 1)
+check_r( 2, Inf, 0, 2)
+
+check_r(-0.2, Inf, 0, 1)
+check_r(-0.2, Inf, 1, 1)
+check_r(-0.2, Inf, 0, 2)
+check_r( 0.0, Inf, 0, 1)
+check_r( 0.0, Inf, 1, 1)
+check_r( 0.0, Inf, 0, 2)
+check_r( 0.2, Inf, 0, 1)
+check_r( 0.2, Inf, 1, 1)
+check_r( 0.2, Inf, 0, 2)
+
+## Right truncation:
+check_r(-Inf, -2, 0, 1)
+check_r(-Inf, -2, 1, 1)
+check_r(-Inf, -2, 0, 2)
+check_r(-Inf, 0, 0, 1)
+check_r(-Inf, 0, 1, 1)
+check_r(-Inf, 0, 0, 2)
+check_r(-Inf, 2, 0, 1)
+check_r(-Inf, 2, 1, 1)
+check_r(-Inf, 2, 0, 2)
+
+check_r(-Inf, -0.2, 0, 1)
+check_r(-Inf, -0.2, 1, 1)
+check_r(-Inf, -0.2, 0, 2)
+check_r(-Inf, 0.0, 0, 1)
+check_r(-Inf, 0.0, 1, 1)
+check_r(-Inf, 0.0, 0, 2)
+check_r(-Inf, 0.2, 0, 1)
+check_r(-Inf, 0.2, 1, 1)
+check_r(-Inf, 0.2, 0, 2)
+
+## Extreme examples:
+check_r(-5, -4, 0, 1)
+
+check_pq <- function(a, b, mean, sd) {
+ for (p in runif(500)) {
+ q <- qtruncnorm(p, a, b, mean, sd)
+ pp <- ptruncnorm(q, a, b, mean, sd)
+ if (abs(p - pp) > 0.00001) {
+ message(sprintf("ptruncnorm(%6.4f, %6.4f, %6.4f, %6.4f, %6.4f) - disagree with qtruncnorm by %f",
+ p, a, b, mean, sd, abs(p - pp)))
+ }
+ }
+}
+
+check_pq(-1, 0, 0, 1)
+check_pq(-1, 1, 0, 1)
+check_pq( 1, 2, 0, 1)
+check_pq(-1, 0, 4, 1)
+check_pq(-1, 1, 4, 1)
+check_pq( 1, 2, 4, 1)
+check_pq(-1, 0, 0, 3)
+check_pq(-1, 1, 0, 3)
+check_pq( 1, 2, 0, 3)
+check_pq(-1, Inf, 0, 1)
+check_pq(-1, Inf, 4, 1)
+check_pq(-1, Inf, 0, 3)
+check_pq(-Inf, 1, 0, 1)
+check_pq(-Inf, 1, 4, 1)
+check_pq(-Inf, 1, 0, 3)
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/r-cran-truncnorm.git
More information about the debian-med-commit
mailing list