The R Project SVN – R – Diff – /trunk/src/nmath/signrank.c



/*
 *  Mathlib : A C Library of Special Functions
 *  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.

        return 1.;
    if (w[0] == 1.)
        return w[k];
 
    w[0] = w[1] = 1.;
    for(j = 2; j < n+1; ++j) {
        int i, end = imin2(j*(j+1)/2, c);
	for(i = end; i >= j; --i)
	    w[i] += w[i-j];
 
    }
 
    return w[k];
}
 

	return(R_D__0);
    x = floor(x + 0.5);
    if ((x < 0) || (x > (n * (n + 1) / 2)))
	return(R_D__0);
 
    int nn = (int) n;
    w_init_maybe(nn);
    d = R_D_exp(log(csignrank((int) x, nn)) - n * M_LN2);
 
    return(d);
}
 
double psignrank(double x, double n, int lower_tail, int log_p)

    if (x < 0.0)
	return(R_DT_0);
    if (x >= n * (n + 1) / 2)
	return(R_DT_1);
 
    int nn = (int) n;
    w_init_maybe(nn);
    f = exp(- n * M_LN2);
    p = 0;
    if (x <= (n * (n + 1) / 4)) {
	for (i = 0; i <= x; i++)
	    p += csignrank(i, nn) * f;
    }
    else {
	x = n * (n + 1) / 2 - x;
	for (i = 0; i < x; i++)
	    p += csignrank(i, nn) * f;
	lower_tail = !lower_tail; /* p = 1 - p; */
    }
 
    return(R_DT_val(p));
} /* psignrank() */
 
double qsignrank(double x, double n, int lower_tail, int log_p)
{
    double f, p;
 
#ifdef IEEE_754
    if (ISNAN(x) || ISNAN(n))
	return(x + n);
#endif

	return(n * (n + 1) / 2);
 
    if(log_p || !lower_tail)
	x = R_DT_qIv(x); /* lower_tail,non-log "p" */
 
    int nn = (int) n;
    w_init_maybe(nn);
    f = exp(- n * M_LN2);
    p = 0;
    int q = 0;
    if (x <= 0.5) {
	x = x - 10 * DBL_EPSILON;
	for (;;) {
	    p += csignrank(q, nn) * f;
	    if (p >= x)
		break;
	    q++;
	}
    }
    else {
	x = 1 - x + 10 * DBL_EPSILON;
	for (;;) {
	    p += csignrank(q, nn) * f;
	    if (p > x) {
		q = (int)(n * (n + 1) / 2 - q);
		break;
	    }
	    q++;
	}
    }

Rev 59167	Rev 59170
Line 1...	Line 1...
1	`/*`	1	`/*`
2	`* Mathlib : A C Library of Special Functions`	2	`* Mathlib : A C Library of Special Functions`
3	`* Copyright (C) 1999-2007 The R Core Team`	3	`* Copyright (C) 1999-2012 The R Core Team`
4	`*`	4	`*`
5	`* This program is free software; you can redistribute it and/or modify`	5	`* This program is free software; you can redistribute it and/or modify`
6	`* it under the terms of the GNU General Public License as published by`	6	`* it under the terms of the GNU General Public License as published by`
7	`* the Free Software Foundation; either version 2 of the License, or`	7	`* the Free Software Foundation; either version 2 of the License, or`
8	`* (at your option) any later version.`	8	`* (at your option) any later version.`
Line 101...	Line 101...
101	`return 1.;`	101	`return 1.;`
102	`if (w[0] == 1.)`	102	`if (w[0] == 1.)`
103	`return w[k];`	103	`return w[k];`
104		104
105	`w[0] = w[1] = 1.;`	105	`w[0] = w[1] = 1.;`
106	`for(j=2; j < n+1; ++j) {`	106	`for(j = 2; j < n+1; ++j) {`
107	`int i, end = imin2(j*(j+1)/2, c);`	107	`int i, end = imin2(j*(j+1)/2, c);`
108	`for(i=end; i >= j; --i) {`	108	`for(i = end; i >= j; --i)`
109	`w[i] += w[i-j];`	109	`w[i] += w[i-j];`
110	`}`	-
111	`}`	110	`}`
112		111
113	`return w[k];`	112	`return w[k];`
114	`}`	113	`}`
115		114
Line 129...	Line 128...
129	`return(R_D__0);`	128	`return(R_D__0);`
130	`x = floor(x + 0.5);`	129	`x = floor(x + 0.5);`
131	`if ((x < 0) \|\| (x > (n * (n + 1) / 2)))`	130	`if ((x < 0) \|\| (x > (n * (n + 1) / 2)))`
132	`return(R_D__0);`	131	`return(R_D__0);`
133		132
-		133	`int nn = (int) n;`
134	`w_init_maybe(n);`	134	`w_init_maybe(nn);`
135	`d = R_D_exp(log(csignrank(x, n)) - n * M_LN2);`	135	`d = R_D_exp(log(csignrank((int) x, nn)) - n * M_LN2);`
136		136
137	`return(d);`	137	`return(d);`
138	`}`	138	`}`
139		139
140	`double psignrank(double x, double n, int lower_tail, int log_p)`	140	`double psignrank(double x, double n, int lower_tail, int log_p)`
Line 154...	Line 154...
154	`if (x < 0.0)`	154	`if (x < 0.0)`
155	`return(R_DT_0);`	155	`return(R_DT_0);`
156	`if (x >= n * (n + 1) / 2)`	156	`if (x >= n * (n + 1) / 2)`
157	`return(R_DT_1);`	157	`return(R_DT_1);`
158		158
-		159	`int nn = (int) n;`
159	`w_init_maybe(n);`	160	`w_init_maybe(nn);`
160	`f = exp(- n * M_LN2);`	161	`f = exp(- n * M_LN2);`
161	`p = 0;`	162	`p = 0;`
162	`if (x <= (n * (n + 1) / 4)) {`	163	`if (x <= (n * (n + 1) / 4)) {`
163	`for (i = 0; i <= x; i++)`	164	`for (i = 0; i <= x; i++)`
164	`p += csignrank(i, n) * f;`	165	`p += csignrank(i, nn) * f;`
165	`}`	166	`}`
166	`else {`	167	`else {`
167	`x = n * (n + 1) / 2 - x;`	168	`x = n * (n + 1) / 2 - x;`
168	`for (i = 0; i < x; i++)`	169	`for (i = 0; i < x; i++)`
169	`p += csignrank(i, n) * f;`	170	`p += csignrank(i, nn) * f;`
170	`lower_tail = !lower_tail; /* p = 1 - p; */`	171	`lower_tail = !lower_tail; /* p = 1 - p; */`
171	`}`	172	`}`
172		173
173	`return(R_DT_val(p));`	174	`return(R_DT_val(p));`
174	`} /* psignrank() */`	175	`} /* psignrank() */`
175		176
176	`double qsignrank(double x, double n, int lower_tail, int log_p)`	177	`double qsignrank(double x, double n, int lower_tail, int log_p)`
177	`{`	178	`{`
178	`double f, p, q;`	179	`double f, p;`
179		180
180	`#ifdef IEEE_754`	181	`#ifdef IEEE_754`
181	`if (ISNAN(x) \|\| ISNAN(n))`	182	`if (ISNAN(x) \|\| ISNAN(n))`
182	`return(x + n);`	183	`return(x + n);`
183	`#endif`	184	`#endif`
Line 195...	Line 196...
195	`return(n * (n + 1) / 2);`	196	`return(n * (n + 1) / 2);`
196		197
197	`if(log_p \|\| !lower_tail)`	198	`if(log_p \|\| !lower_tail)`
198	`x = R_DT_qIv(x); /* lower_tail,non-log "p" */`	199	`x = R_DT_qIv(x); /* lower_tail,non-log "p" */`
199		200
-		201	`int nn = (int) n;`
200	`w_init_maybe(n);`	202	`w_init_maybe(nn);`
201	`f = exp(- n * M_LN2);`	203	`f = exp(- n * M_LN2);`
202	`p = 0;`	204	`p = 0;`
203	`q = 0;`	205	`int q = 0;`
204	`if (x <= 0.5) {`	206	`if (x <= 0.5) {`
205	`x = x - 10 * DBL_EPSILON;`	207	`x = x - 10 * DBL_EPSILON;`
206	`for (;;) {`	208	`for (;;) {`
207	`p += csignrank(q, n) * f;`	209	`p += csignrank(q, nn) * f;`
208	`if (p >= x)`	210	`if (p >= x)`
209	`break;`	211	`break;`
210	`q++;`	212	`q++;`
211	`}`	213	`}`
212	`}`	214	`}`
213	`else {`	215	`else {`
214	`x = 1 - x + 10 * DBL_EPSILON;`	216	`x = 1 - x + 10 * DBL_EPSILON;`
215	`for (;;) {`	217	`for (;;) {`
216	`p += csignrank(q, n) * f;`	218	`p += csignrank(q, nn) * f;`
217	`if (p > x) {`	219	`if (p > x) {`
218	`q = n * (n + 1) / 2 - q;`	220	`q = (int)(n * (n + 1) / 2 - q);`
219	`break;`	221	`break;`
220	`}`	222	`}`
221	`q++;`	223	`q++;`
222	`}`	224	`}`
223	`}`	225	`}`

The R Project SVN R

(root)/trunk/src/nmath/signrank.c – Rev 59167 → 59170