1 /*
2 * R : A Computer Language for Statistical Data Analysis
3 * Copyright (C) 1999-2016 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 * https://www.R-project.org/Licenses/
18 *
19 */
20
21 /* ansari.c
22 Compute the exact distribution of the Ansari-Bradley test statistic.
23 */
24
25 #include <string.h>
26 #include <R.h>
27 #include <math.h> // for floor
28 #include <Rmath.h> /* uses choose() */
29 #include "stats.h"
30
31 static double ***
w_init(int m,int n)32 w_init(int m, int n)
33 {
34 int i;
35 double ***w;
36
37 w = (double ***) R_alloc(m + 1, sizeof(double **));
38 memset(w, '\0', (m+1) * sizeof(double**));
39 for (i = 0; i <= m; i++) {
40 w[i] = (double**) R_alloc(n + 1, sizeof(double *));
41 memset(w[i], '\0', (n+1) * sizeof(double*));
42 }
43 return(w);
44 }
45
46
47 static double
cansari(int k,int m,int n,double *** w)48 cansari(int k, int m, int n, double ***w)
49 {
50 int i, l, u;
51
52 l = (m + 1) * (m + 1) / 4;
53 u = l + m * n / 2;
54
55 if ((k < l) || (k > u))
56 return(0);
57
58 if (w[m][n] == 0) {
59 w[m][n] = (double *) R_alloc(u + 1, sizeof(double));
60 memset(w[m][n], '\0', (u + 1) * sizeof(double));
61 for (i = 0; i <= u; i++)
62 w[m][n][i] = -1;
63 }
64
65 if (w[m][n][k] < 0) {
66 if (m == 0)
67 w[m][n][k] = (k == 0);
68 else if (n == 0)
69 w[m][n][k] = (k == l);
70 else
71 w[m][n][k] = cansari(k, m, n - 1, w)
72 + cansari(k - (m + n + 1) / 2, m - 1, n, w);
73 }
74
75 return(w[m][n][k]);
76 }
77
78
79 static void
pansari(int len,double * Q,double * P,int m,int n)80 pansari(int len, double *Q, double *P, int m, int n)
81 {
82 int i, j, l, u;
83 double c, p, q;
84 double ***w;
85
86 w = w_init(m, n);
87 l = (m + 1) * (m + 1) / 4;
88 u = l + m * n / 2;
89 c = choose(m + n, m);
90 for (i = 0; i < len; i++) {
91 q = floor(Q[i] + 1e-7);
92 if (q < l)
93 P[i] = 0;
94 else if (q > u)
95 P[i] = 1;
96 else {
97 p = 0;
98 for (j = l; j <= q; j++) p += cansari(j, m, n, w);
99 P[i] = p / c;
100 }
101 }
102 }
103
104 static void
qansari(int len,double * P,double * Q,int m,int n)105 qansari(int len, double *P, double *Q, int m, int n)
106 {
107 int i, l, u;
108 double c, p, xi;
109 double ***w;
110
111 w = w_init(m, n);
112 l = (m + 1) * (m + 1) / 4;
113 u = l + m * n / 2;
114 c = choose(m + n, m);
115 for (i = 0; i < len; i++) {
116 xi = P[i];
117 if(xi < 0 || xi > 1)
118 error(_("probabilities outside [0,1] in qansari()"));
119 if(xi == 0)
120 Q[i] = l;
121 else if(xi == 1)
122 Q[i] = u;
123 else {
124 p = 0.;
125 int q = 0;
126 for(;;) {
127 p += cansari(q, m, n, w) / c;
128 if (p >= xi) break;
129 q++;
130 }
131 Q[i] = q;
132 }
133 }
134 }
135
136 #include <Rinternals.h>
pAnsari(SEXP q,SEXP sm,SEXP sn)137 SEXP pAnsari(SEXP q, SEXP sm, SEXP sn)
138 {
139 int m = asInteger(sm), n = asInteger(sn);
140 q = PROTECT(coerceVector(q, REALSXP));
141 int len = LENGTH(q);
142 SEXP p = PROTECT(allocVector(REALSXP, len));
143 pansari(len, REAL(q), REAL(p), m, n);
144 UNPROTECT(2);
145 return p;
146 }
147
qAnsari(SEXP p,SEXP sm,SEXP sn)148 SEXP qAnsari(SEXP p, SEXP sm, SEXP sn)
149 {
150 int m = asInteger(sm), n = asInteger(sn);
151 p = PROTECT(coerceVector(p, REALSXP));
152 int len = LENGTH(p);
153 SEXP q = PROTECT(allocVector(REALSXP, len));
154 qansari(len, REAL(p), REAL(q), m, n);
155 UNPROTECT(2);
156 return q;
157 }
158