[Distr-commits] r907 - in branches/distr-2.5/pkg/distr: . inst src

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Wed Sep 11 18:47:08 CEST 2013


Author: ruckdeschel
Date: 2013-09-11 18:47:07 +0200 (Wed, 11 Sep 2013)
New Revision: 907

Added:
   branches/distr-2.5/pkg/distr/src/
   branches/distr-2.5/pkg/distr/src/ks.c
Modified:
   branches/distr-2.5/pkg/distr/NAMESPACE
   branches/distr-2.5/pkg/distr/inst/NEWS
Log:
[distr] backported changes in trunc to devel

Modified: branches/distr-2.5/pkg/distr/NAMESPACE
===================================================================
--- branches/distr-2.5/pkg/distr/NAMESPACE	2013-09-11 16:32:15 UTC (rev 906)
+++ branches/distr-2.5/pkg/distr/NAMESPACE	2013-09-11 16:47:07 UTC (rev 907)
@@ -1,3 +1,4 @@
+useDynLib("distr")
 import("methods")
 import("stats")
 importFrom("graphics", "plot")

Modified: branches/distr-2.5/pkg/distr/inst/NEWS
===================================================================
--- branches/distr-2.5/pkg/distr/inst/NEWS	2013-09-11 16:32:15 UTC (rev 906)
+++ branches/distr-2.5/pkg/distr/inst/NEWS	2013-09-11 16:47:07 UTC (rev 907)
@@ -20,7 +20,11 @@
 + added .Rbuildignore
 + exported some routines which had been internal so far to
   avoid calls by :::
-
++ in order to avoid calls to non-exported C-Code from pkg stats
+  (for Kolmogorov-distribution) we integrated R-Core's C-code from  
+            /src/library/stats/src/ks.c rev60573
+   (with small changes to make it self-contained) to src folder of this pkg;
+    in addition, we now use the shakier .C - interface again...
 BUGFIXES:
 + fixed .makedotsPt -issue discovered by Gerald --> thx
 + moved generics to "distribution" and "samplesize" to pkg distr to avoid 

Added: branches/distr-2.5/pkg/distr/src/ks.c
===================================================================
--- branches/distr-2.5/pkg/distr/src/ks.c	                        (rev 0)
+++ branches/distr-2.5/pkg/distr/src/ks.c	2013-09-11 16:47:07 UTC (rev 907)
@@ -0,0 +1,261 @@
+/*
+ * taken from R Core: /src/library/stats/src/ks.c rev60573
+ */  
+/*
+ *  R : A Computer Language for Statistical Data Analysis
+ *  Copyright (C) 1999-2012   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/
+ */
+
+/* ks.c
+   Compute the asymptotic distribution of the one- and two-sample
+   two-sided Kolmogorov-Smirnov statistics, and the exact distributions
+   in the two-sided one-sample and two-sample cases.
+*/
+
+#include <stdlib.h>
+#include <math.h>
+#define PI 3.141592653589793
+#ifndef M_1_SQRT_2PI
+#define M_1_SQRT_2PI	0.398942280401432677939946059934	/* 1/sqrt(2pi) */
+#endif
+#ifndef M_PI_2
+#define M_PI_2		1.570796326794896619231321691640	/* pi/2 */
+#endif
+
+#ifndef M_PI_4
+#define M_PI_4		0.785398163397448309615660845820	/* pi/4 */
+#endif
+
+/*
+ * taken from R Core: /src/include/R_ext/RS.h rev60573
+ */  
+extern void *R_chk_calloc(size_t, size_t);
+extern void *R_chk_realloc(void *, size_t);
+extern void R_chk_free(void *);
+
+/*
+ * taken from R Core: /src/main/memory.c rev60573
+ */  
+void *R_chk_calloc(size_t nelem, size_t elsize)
+{
+    void *p;
+    p = calloc(nelem, elsize);
+    return(p);
+}
+
+void *R_chk_realloc(void *ptr, size_t size)
+{
+    void *p;
+    /* Protect against broken realloc */
+    if(ptr) p = realloc(ptr, size); else p = malloc(size);
+    return(p);
+}
+
+void R_chk_free(void *ptr)
+{
+    /* S-PLUS warns here, but there seems no reason to do so */
+    /* if(!ptr) warning("attempt to free NULL pointer by Free"); */
+    if(ptr) free(ptr); /* ANSI C says free has no effect on NULL, but
+			  better to be safe here */
+}
+
+#define Calloc(n, t)   (t *) R_chk_calloc( (size_t) (n), sizeof(t) )
+#define Realloc(p,n,t) (t *) R_chk_realloc( (void *)(p), (size_t)((n) * sizeof(t)) )
+#define Free(p)        (R_chk_free( (void *)(p) ), (p) = NULL)
+
+/*static void pkolmogorov2x(double *x, int *n);
+static void pkstwo(int n, double *x, double tol);
+*/
+static double K(int n, double d);
+static void m_multiply(double *A, double *B, double *C, int m);
+static void m_power(double *A, int eA, double *V, int *eV, int m, int n);
+
+/* Two-sample two-sided asymptotic distribution */
+void
+pkstwo(int *n, double *x, double *tol)
+{
+/* x[1:n] is input and output
+ *
+ * Compute
+ *   \sum_{k=-\infty}^\infty (-1)^k e^{-2 k^2 x^2}
+ *   = 1 + 2 \sum_{k=1}^\infty (-1)^k e^{-2 k^2 x^2}
+ *   = \frac{\sqrt{2\pi}}{x} \sum_{k=1}^\infty \exp(-(2k-1)^2\pi^2/(8x^2))
+ *
+ * See e.g. J. Durbin (1973), Distribution Theory for Tests Based on the
+ * Sample Distribution Function.  SIAM.
+ *
+ * The 'standard' series expansion obviously cannot be used close to 0;
+ * we use the alternative series for x < 1, and a rather crude estimate
+ * of the series remainder term in this case, in particular using that
+ * ue^(-lu^2) \le e^(-lu^2 + u) \le e^(-(l-1)u^2 - u^2+u) \le e^(-(l-1))
+ * provided that u and l are >= 1.
+ *
+ * (But note that for reasonable tolerances, one could simply take 0 as
+ * the value for x < 0.2, and use the standard expansion otherwise.)
+ *
+ */
+    double new, old, s, w, z;
+    int i, k, k_max;
+
+    k_max = (int) sqrt(2 - log(*tol));
+
+    for(i = 0; i < *n; i++) {
+	if(x[i] < 1) {
+	    z = - (M_PI_2 * M_PI_4) / (x[i] * x[i]);
+	    w = log(x[i]);
+	    s = 0;
+	    for(k = 1; k < k_max; k += 2) {
+		s += exp(k * k * z - w);
+	    }
+	    x[i] = s / M_1_SQRT_2PI;
+	}
+	else {
+	    z = -2 * x[i] * x[i];
+	    s = -1;
+	    k = 1;
+	    old = 0;
+	    new = 1;
+	    while(fabs(old - new) > *tol) {
+		old = new;
+		new += 2 * s * exp(z * k * k);
+		s *= -1;
+		k++;
+	    }
+	    x[i] = new;
+	}
+    }
+}
+
+static double
+K(int n, double d)
+{
+    /* Compute Kolmogorov's distribution.
+       Code published in
+	 George Marsaglia and Wai Wan Tsang and Jingbo Wang (2003),
+	 "Evaluating Kolmogorov's distribution".
+	 Journal of Statistical Software, Volume 8, 2003, Issue 18.
+	 URL: http://www.jstatsoft.org/v08/i18/.
+    */
+
+   int k, m, i, j, g, eH, eQ;
+   double h, s, *H, *Q;
+
+   /* 
+      The faster right-tail approximation is omitted here.
+      s = d*d*n; 
+      if(s > 7.24 || (s > 3.76 && n > 99)) 
+          return 1-2*exp(-(2.000071+.331/sqrt(n)+1.409/n)*s);
+   */
+   k = (int) (n * d) + 1;
+   m = 2 * k - 1;
+   h = k - n * d;
+   H = (double*) Calloc(m * m, double);
+   Q = (double*) Calloc(m * m, double);
+   for(i = 0; i < m; i++)
+       for(j = 0; j < m; j++)
+	   if(i - j + 1 < 0)
+	       H[i * m + j] = 0;
+	   else
+	       H[i * m + j] = 1;
+   for(i = 0; i < m; i++) {
+       H[i * m] -= pow(h, i + 1);
+       H[(m - 1) * m + i] -= pow(h, (m - i));
+   }
+   H[(m - 1) * m] += ((2 * h - 1 > 0) ? pow(2 * h - 1, m) : 0);
+   for(i = 0; i < m; i++)
+       for(j=0; j < m; j++)
+	   if(i - j + 1 > 0)
+	       for(g = 1; g <= i - j + 1; g++)
+		   H[i * m + j] /= g;
+   eH = 0;
+   m_power(H, eH, Q, &eQ, m, n);
+   s = Q[(k - 1) * m + k - 1];
+   for(i = 1; i <= n; i++) {
+       s = s * i / n;
+       if(s < 1e-140) {
+	   s *= 1e140;
+	   eQ -= 140;
+       }
+   }
+   s *= pow(10., eQ);
+   Free(H);
+   Free(Q);
+   return(s);
+}
+
+static void
+m_multiply(double *A, double *B, double *C, int m)
+{
+    /* Auxiliary routine used by K().
+       Matrix multiplication.
+    */
+    int i, j, k;
+    double s;
+    for(i = 0; i < m; i++)
+	for(j = 0; j < m; j++) {
+	    s = 0.;
+	    for(k = 0; k < m; k++)
+		s+= A[i * m + k] * B[k * m + j];
+	    C[i * m + j] = s;
+	}
+}
+
+static void
+m_power(double *A, int eA, double *V, int *eV, int m, int n)
+{
+    /* Auxiliary routine used by K().
+       Matrix power.
+    */
+    double *B;
+    int eB , i;
+
+    if(n == 1) {
+	for(i = 0; i < m * m; i++)
+	    V[i] = A[i];
+	*eV = eA;
+	return;
+    }
+    m_power(A, eA, V, eV, m, n / 2);
+    B = (double*) Calloc(m * m, double);
+    m_multiply(V, V, B, m);
+    eB = 2 * (*eV);
+    if((n % 2) == 0) {
+	for(i = 0; i < m * m; i++)
+	    V[i] = B[i];
+	*eV = eB;
+    }
+    else {
+	m_multiply(A, B, V, m);
+	*eV = eA + eB;
+    }
+    if(V[(m / 2) * m + (m / 2)] > 1e140) {
+	for(i = 0; i < m * m; i++)
+	    V[i] = V[i] * 1e-140;
+	*eV += 140;
+    }
+    Free(B);
+}
+/*
+ * taken from R Core: /src/library/stats/src/ks.c rev56090
+*/
+void pkolmogorov2x(double *x, int *n)
+{
+    /* x is input and output. */
+
+    *x = K(*n, *x);
+}
+



More information about the Distr-commits mailing list