Commit 09af4c0f authored by Reinhard Furrer's avatar Reinhard Furrer

update with c function

parent 90afff4e
......@@ -2,3 +2,8 @@
.Rhistory
.RData
.Ruserdata
src/*.o
src/*.so
*.*~
......@@ -9,3 +9,6 @@ importFrom('wmtsa', 'wavCWT')
export("prsim")
export("prsim.wave")
useDynLib(PRSim, .registration = TRUE)
......@@ -25,7 +25,7 @@ ks_test <- function (x, y)
p[is.na(x)] <- NA
IND <- which(!is.na(x) & (x > 0))
if (length(IND))
p[IND] <- .Call(stats:::C_pKS2, p = x[IND], tol)
p[IND] <- .Call(pks2, p = x[IND], tol, PACKAGE="PRSim")
p
}
PVAL <- 1 - pkstwo(sqrt(n) * STATISTIC)
......
/* Reinhard Furrer, fall 2019, based on spam's version and 'Writing R extensions'. */
#include <R_ext/RS.h>
#include <stdlib.h> // for NULL
#include <R_ext/Rdynload.h>
extern void pks2 ( void *, void *);
/* to get all functions:
nm -g ./PRSim.so | grep " T "
*/
/* .Call calls */
static const R_CallMethodDef CallEntries[] = {
{"pks2", (DL_FUNC) &pks2, 2},
{NULL, NULL, 0}
};
void R_init_PRSim(DllInfo *dll)
{
R_registerRoutines(dll, NULL, CallEntries, NULL, NULL);
R_useDynamicSymbols(dll, FALSE);
}
/*
* R : A Computer Language for Statistical Data Analysis
* Copyright (C) 1999-2016 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
* https://www.R-project.org/Licenses/
*/
/* this is a stripped down version of ks.c
r-source/blob/trunk/src/library/stats/src/ks.c
Compute the asymptotic distribution of the two-sample
two-sided Kolmogorov-Smirnov statistics.
*/
#include <math.h>
#include <R.h>
#include <Rinternals.h>
#include <Rmath.h> /* constants */
/* Two-sample two-sided asymptotic distribution */
static 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;
}
}
}
/* Two-sample two-sided asymptotic distribution */
SEXP pks2(SEXP statistic, SEXP stol)
{
int n = LENGTH(statistic);
double tol = asReal(stol);
SEXP ans = duplicate(statistic);
pkstwo(n, REAL(ans), tol);
return ans;
}
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment