[Desire-commits] r17 - in packages/truncnorm: . src

[noreply at r-forge.r-project.org]

Author: olafm
Date: 2008-05-21 10:12:00 +0200 (Wed, 21 May 2008)
New Revision: 17

Modified:
   packages/truncnorm/DESCRIPTION
   packages/truncnorm/src/truncnorm.c
Log:
* FIX: Misuse of pnorm()
* Rework constant calculation


Modified: packages/truncnorm/DESCRIPTION
===================================================================
--- packages/truncnorm/DESCRIPTION	2008-05-20 12:28:40 UTC (rev 16)
+++ packages/truncnorm/DESCRIPTION	2008-05-21 08:12:00 UTC (rev 17)
@@ -1,6 +1,6 @@
 Package: truncnorm
-Version: 0.9.2
-Date: 2008-05-19
+Version: 0.9.3
+Date: 2008-05-20
 Title: Trunctated Normal Distribution
 Author:
   Heike Trautmann <trautmann at statistik.uni-dortmund.de>

Modified: packages/truncnorm/src/truncnorm.c
===================================================================
--- packages/truncnorm/src/truncnorm.c	2008-05-20 12:28:40 UTC (rev 16)
+++ packages/truncnorm/src/truncnorm.c	2008-05-21 08:12:00 UTC (rev 17)
@@ -38,7 +38,7 @@
 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;
+  double *x, a, b, mean, sd, *ret;
 
   UNPACK_REAL_VECTOR(args, s_x, x);
   UNPACK_REAL_VALUE(args, s_a, a);
@@ -48,13 +48,14 @@
 
   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);
+  const double c1 = pnorm((a-mean)/sd, 0.0, 1.0, TRUE, FALSE);
+  const double c2 = pnorm((b-mean)/sd, 0.0, 1.0, TRUE, FALSE);
+  const double c3 = sd * (c2 - c1);
   for (i = 0; i < n; ++i) {
-    const double t = x[i];
-    if (a <= t && t <= b) { /* In range: */
-      c1 = dnorm(t, mean, sd, FALSE);
-      ret[i] = - c1 / c2;
+    const double cx = x[i];
+    if (a <= cx && cx <= b) { /* In range: */
+      const double c4 = dnorm((cx - mean)/sd, 0.0, 1.0, FALSE);
+      ret[i] = c4 / c3;
     } else { /* Truncated: */
       ret[i] = 0.0;
     }
@@ -78,8 +79,8 @@
   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);
+  ca = pnorm(a, mean, sd, TRUE, FALSE);
+  cb = pnorm(b, mean, sd, TRUE, FALSE);
   
   for (i = 0; i < n; ++i) {
     if (q[i] < a) {
@@ -87,7 +88,7 @@
     } else if (q[i] > b) {
       ret[i] = 1.0;
     } else {
-      const double cx = pnorm(q[i], mean, sd, FALSE, FALSE);
+      const double cx = pnorm(q[i], mean, sd, TRUE, FALSE);
       ret[i] = (cx - ca) / (cb - ca);
     }
   }
@@ -146,8 +147,8 @@
   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);
+    const double Ca = pnorm(a[i], mean[i], sd[i], TRUE, FALSE);
+    const double Cb = pnorm(b[i], mean[i], sd[i], TRUE, FALSE);
     
     ret[i] = mean[i] + sd[i] * ((ca - cb) / (Cb - Ca));
   } 
@@ -180,8 +181,8 @@
     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 Ca = pnorm(am, 0.0, 1.0, TRUE, FALSE);
+    const double Cb = pnorm(bm, 0.0, 1.0, TRUE, FALSE);
     const double v  = sd[i] * sd[i];
     
     const double d = Cb - Ca;