1 /*
2 * Mathlib : A C Library of Special Functions
3 * Copyright (C) 1999-2012 The R Core Team
4 *
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
7 * the Free Software Foundation; either version 2 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, a copy is available at
17 * http://www.r-project.org/Licenses/
18 *
19 * SYNOPSIS
20 *
21 * #include <Rmath.h>
22 * double dsignrank(double x, double n, int give_log)
23 * double psignrank(double x, double n, int lower_tail, int log_p)
24 * double qsignrank(double x, double n, int lower_tail, int log_p)
25 * double rsignrank(double n)
26 *
27 * DESCRIPTION
28 *
29 * dsignrank The density of the Wilcoxon Signed Rank distribution.
30 * psignrank The distribution function of the Wilcoxon Signed Rank
31 * distribution.
32 * qsignrank The quantile function of the Wilcoxon Signed Rank
33 * distribution.
34 * rsignrank Random variates from the Wilcoxon Signed Rank
35 * distribution.
36 */
37
38 #include "nmath.h"
39 #include "dpq.h"
40
41 static double *w;
42 static int allocated_n;
43
44 static void
w_free(void)45 w_free(void)
46 {
47 if (!w) return;
48
49 free((void *) w);
50 w = 0;
51 allocated_n = 0;
52 }
53
signrank_free(void)54 void signrank_free(void)
55 {
56 w_free();
57 }
58
59 static void
w_init_maybe(int n)60 w_init_maybe(int n)
61 {
62 int u, c;
63
64 u = n * (n + 1) / 2;
65 c = (u / 2);
66
67 if (w) {
68 if(n != allocated_n) {
69 w_free();
70 }
71 else return;
72 }
73
74 if(!w) {
75 w = (double *) calloc((size_t) c + 1, sizeof(double));
76 #ifdef MATHLIB_STANDALONE
77 if (!w) MATHLIB_ERROR("%s", _("signrank allocation error"));
78 #endif
79 allocated_n = n;
80 }
81 }
82
83 static double
csignrank(int k,int n)84 csignrank(int k, int n)
85 {
86 int c, u, j;
87
88 #ifndef MATHLIB_STANDALONE
89 R_CheckUserInterrupt();
90 #endif
91
92 u = n * (n + 1) / 2;
93 c = (u / 2);
94
95 if (k < 0 || k > u)
96 return 0;
97 if (k > c)
98 k = u - k;
99
100 if (n == 1)
101 return 1.;
102 if (w[0] == 1.)
103 return w[k];
104
105 w[0] = w[1] = 1.;
106 for(j = 2; j < n+1; ++j) {
107 int i, end = imin2(j*(j+1)/2, c);
108 for(i = end; i >= j; --i)
109 w[i] += w[i-j];
110 }
111
112 return w[k];
113 }
114
dsignrank(double x,double n,int give_log)115 double dsignrank(double x, double n, int give_log)
116 {
117 double d;
118
119 #ifdef IEEE_754
120 /* NaNs propagated correctly */
121 if (ISNAN(x) || ISNAN(n)) return(x + n);
122 #endif
123 n = floor(n + 0.5);
124 if (n <= 0)
125 ML_ERR_return_NAN;
126
127 if (fabs(x - floor(x + 0.5)) > 1e-7)
128 return(R_D__0);
129 x = floor(x + 0.5);
130 if ((x < 0) || (x > (n * (n + 1) / 2)))
131 return(R_D__0);
132
133 int nn = (int) n;
134 w_init_maybe(nn);
135 d = R_D_exp(log(csignrank((int) x, nn)) - n * M_LN2);
136
137 return(d);
138 }
139
psignrank(double x,double n,int lower_tail,int log_p)140 double psignrank(double x, double n, int lower_tail, int log_p)
141 {
142 int i;
143 double f, p;
144
145 #ifdef IEEE_754
146 if (ISNAN(x) || ISNAN(n))
147 return(x + n);
148 #endif
149 if (!R_FINITE(n)) ML_ERR_return_NAN;
150 n = floor(n + 0.5);
151 if (n <= 0) ML_ERR_return_NAN;
152
153 x = floor(x + 1e-7);
154 if (x < 0.0)
155 return(R_DT_0);
156 if (x >= n * (n + 1) / 2)
157 return(R_DT_1);
158
159 int nn = (int) n;
160 w_init_maybe(nn);
161 f = exp(- n * M_LN2);
162 p = 0;
163 if (x <= (n * (n + 1) / 4)) {
164 for (i = 0; i <= x; i++)
165 p += csignrank(i, nn) * f;
166 }
167 else {
168 x = n * (n + 1) / 2 - x;
169 for (i = 0; i < x; i++)
170 p += csignrank(i, nn) * f;
171 lower_tail = !lower_tail; /* p = 1 - p; */
172 }
173
174 return(R_DT_val(p));
175 } /* psignrank() */
176
qsignrank(double x,double n,int lower_tail,int log_p)177 double qsignrank(double x, double n, int lower_tail, int log_p)
178 {
179 double f, p;
180
181 #ifdef IEEE_754
182 if (ISNAN(x) || ISNAN(n))
183 return(x + n);
184 #endif
185 if (!R_FINITE(x) || !R_FINITE(n))
186 ML_ERR_return_NAN;
187 R_Q_P01_check(x);
188
189 n = floor(n + 0.5);
190 if (n <= 0)
191 ML_ERR_return_NAN;
192
193 if (x == R_DT_0)
194 return(0);
195 if (x == R_DT_1)
196 return(n * (n + 1) / 2);
197
198 if(log_p || !lower_tail)
199 x = R_DT_qIv(x); /* lower_tail,non-log "p" */
200
201 int nn = (int) n;
202 w_init_maybe(nn);
203 f = exp(- n * M_LN2);
204 p = 0;
205 int q = 0;
206 if (x <= 0.5) {
207 x = x - 10 * DBL_EPSILON;
208 for (;;) {
209 p += csignrank(q, nn) * f;
210 if (p >= x)
211 break;
212 q++;
213 }
214 }
215 else {
216 x = 1 - x + 10 * DBL_EPSILON;
217 for (;;) {
218 p += csignrank(q, nn) * f;
219 if (p > x) {
220 q = (int)(n * (n + 1) / 2 - q);
221 break;
222 }
223 q++;
224 }
225 }
226
227 return(q);
228 }
229
rsignrank(double n)230 double rsignrank(double n)
231 {
232 int i, k;
233 double r;
234
235 #ifdef IEEE_754
236 /* NaNs propagated correctly */
237 if (ISNAN(n)) return(n);
238 #endif
239 n = floor(n + 0.5);
240 if (n < 0) ML_ERR_return_NAN;
241
242 if (n == 0)
243 return(0);
244 r = 0.0;
245 k = (int) n;
246 for (i = 0; i < k; ) {
247 r += (++i) * floor(unif_rand() + 0.5);
248 }
249 return(r);
250 }
251