I took some time and liberty and tried to improve
existing implementation of SignRank functions
in R. (dsignrank, ...)
As I have seen they've been based on csignrank.
So I modified csignrank and, I believe,
improved calculation time and memory efficiency.
The idea is basically the same. I use the same recursion
as original author used with one slight modification.
I am generating Wilcoxon SignRank density from the
beginning (for n=2,3,...) with the help of recursion formula.
What is changed?
There is no need for SINGRANK_MAX in src/nmath/nmath.h anymore.
Only functions:
static void w_free()
void signrank_free()
static void w_init_maybe(int n)
static double csignrank(int k, int n)
in src/nmath/signrank.c have been altered.
There was no change to dsignrank, psignrank, ...
I've tried to make as little changes as possible.
So, to compile with new functions only src/nmath/signrank.c
needs to be changed with the one given in attachment.
I hope it really is an improvement.
With respect,
--
Ivo Ugrina
ICQ: 47508335 | www.iugrina.com
---
baza matematickih pojmova
http://baza.iugrina.com
---
anime, manga, Japan fanzin
http://yoshi.iugrina.com
/*
* Mathlib : A C Library of Special Functions
* Copyright (C) 1999-2007 The R Development 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/
*
* SYNOPSIS
*
*#include Rmath.h
*double dsignrank(double x, double n, int give_log)
*double psignrank(double x, double n, int lower_tail, int log_p)
*double qsignrank(double x, double n, int lower_tail, int log_p)
*double rsignrank(double n)
*
* DESCRIPTION
*
*dsignrankThe density of the Wilcoxon Signed Rank distribution.
*psignrankThe distribution function of the Wilcoxon Signed Rank
* distribution.
*qsignrankThe quantile function of the Wilcoxon Signed Rank
* distribution.
*rsignrankRandom variates from the Wilcoxon Signed Rank
* distribution.
*/
#include nmath.h
#include dpq.h
static double *w;
static int allocated_n;
static void
w_free()
{
if( !w) return;
free( (void *) w);
w = 0;
allocated_n = 0;
}
void signrank_free()
{
w_free();
}
static void
w_init_maybe(int n)
{
int u, c;
u = n * (n + 1) / 2;
c = (int) (u / 2);
if( w){
if( n != allocated_n){
w_free();
}
else return;
}
if( !w){
w = (double *) calloc(c + 1, sizeof(double));
#ifdef MATHLIB_STANDALONE
if (!w) MATHLIB_ERROR(%s, _(signrank allocation error));
#endif
allocated_n = n;
}
}
static double
csignrank(int k, int n)
{
int c, u, i, j, end;
#ifndef MATHLIB_STANDALONE
R_CheckUserInterrupt();
#endif
u = n * (n + 1) / 2;
c = (int) (u / 2);
if ((k 0) || (k u))
return 0;
if (k c)
k = u - k;
if ( n == 1)
return 1.0;
if ( w[0]==1.0 )
return w[k];
w[0]=w[1]=1.0;
for( j=2; jn+1; ++j){
end=j*(j+1)/2;
end= ( endc )?end:c;
for( i=end; i=j; --i){
w[i] += w[i-j];
}
}
return w[k];
}
double dsignrank(double x, double n, int give_log)
{
double d;
#ifdef IEEE_754
/* NaNs propagated correctly */
if (ISNAN(x) || ISNAN(n)) return(x + n);
#endif
n = floor(n + 0.5);
if (n = 0)
ML_ERR_return_NAN;
if (fabs(x - floor(x + 0.5)) 1e-7)
return(R_D__0);
x = floor(x + 0.5);
if ((x 0) || (x (n * (n + 1) / 2)))
return(R_D__0);
w_init_maybe(n);
d = R_D_exp(log(csignrank(x, n)) - n * M_LN2);
return(d);
}
double psignrank(double x, double n, int lower_tail, int log_p)
{
int i;
double f, p;
#ifdef IEEE_754
if (ISNAN(x) || ISNAN(n))
return(x + n);
#endif
if (!R_FINITE(n)) ML_ERR_return_NAN;
n = floor(n + 0.5);
if (n = 0) ML_ERR_return_NAN;
x = floor(x + 1e-7);
if (x 0.0)
return(R_DT_0);
if (x = n * (n + 1) / 2)
return(R_DT_1);
w_init_maybe(n);
f = exp(- n * M_LN2);
p = 0;
if (x = (n * (n + 1) / 4)) {
for (i = 0; i = x; i++)
p += csignrank(i, n) * f;
}
else {
x = n * (n + 1) / 2 - x;
for (i = 0; i x; i++)
p += csignrank(i, n) * f;
lower_tail
