1 /* PSPP - a program for statistical analysis.
2 Copyright (C) 2011 Free Software Foundation, Inc.
3
4 This program is free software: you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation, either version 3 of the License, or
7 (at your option) any later version.
8
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with this program. If not, see <http://www.gnu.org/licenses/>.
16 */
17
18 /* This file is taken from the R project source code, and modified.
19 The original copyright notice is reproduced below: */
20
21 /*
22 * Mathlib : A C Library of Special Functions
23 * Copyright (C) 1998 Ross Ihaka
24 * Copyright (C) 2000--2005 The R Development Core Team
25 * based in part on AS70 (C) 1974 Royal Statistical Society
26 *
27 * This program is free software; you can redistribute it and/or modify
28 * it under the terms of the GNU General Public License as published by
29 * the Free Software Foundation; either version 2 of the License, or
30 * (at your option) any later version.
31 *
32 * This program is distributed in the hope that it will be useful,
33 * but WITHOUT ANY WARRANTY; without even the implied warranty of
34 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
35 * GNU General Public License for more details.
36 *
37 * You should have received a copy of the GNU General Public License
38 * along with this program; if not, a copy is available at
39 * http://www.r-project.org/Licenses/
40 *
41 * SYNOPSIS
42 *
43 * #include <Rmath.h>
44 * double qtukey(p, rr, cc, df, lower_tail, log_p);
45 *
46 * DESCRIPTION
47 *
48 * Computes the quantiles of the maximum of rr studentized
49 * ranges, each based on cc means and with df degrees of freedom
50 * for the standard error, is less than q.
51 *
52 * The algorithm is based on that of the reference.
53 *
54 * REFERENCE
55 *
56 * Copenhaver, Margaret Diponzio & Holland, Burt S.
57 * Multiple comparisons of simple effects in
58 * the two-way analysis of variance with fixed effects.
59 * Journal of Statistical Computation and Simulation,
60 * Vol.30, pp.1-15, 1988.
61 */
62
63 #include <config.h>
64
65 #include "tukey.h"
66
67 #include <assert.h>
68 #include <math.h>
69
70 #define TRUE (1)
71 #define FALSE (0)
72
73 #define ML_POSINF (1.0 / 0.0)
74 #define ML_NEGINF (-1.0 / 0.0)
75
76 #define R_D_Lval(p) (lower_tail ? (p) : (0.5 - (p) + 0.5)) /* p */
77
78 #define R_DT_qIv(p) (log_p ? (lower_tail ? exp(p) : - expm1(p)) \
79 : R_D_Lval(p))
80
81
fmax2(double x,double y)82 static double fmax2(double x, double y)
83 {
84 if (isnan(x) || isnan(y))
85 return x + y;
86 return (x < y) ? y : x;
87 }
88
89
90 #define R_Q_P01_boundaries(p, _LEFT_, _RIGHT_) \
91 if (log_p) { \
92 assert (p <= 0); \
93 if(p == 0) /* upper bound*/ \
94 return lower_tail ? _RIGHT_ : _LEFT_; \
95 if(p == ML_NEGINF) \
96 return lower_tail ? _LEFT_ : _RIGHT_; \
97 } \
98 else { /* !log_p */ \
99 assert (p >= 0 && p <= 1); \
100 if(p == 0) \
101 return lower_tail ? _LEFT_ : _RIGHT_; \
102 if(p == 1) \
103 return lower_tail ? _RIGHT_ : _LEFT_; \
104 }
105
106
107 /* qinv() :
108 * this function finds percentage point of the studentized range
109 * which is used as initial estimate for the secant method.
110 * function is adapted from portion of algorithm as 70
111 * from applied statistics (1974) ,vol. 23, no. 1
112 * by odeh, r. e. and evans, j. o.
113 *
114 * p = percentage point
115 * c = no. of columns or treatments
116 * v = degrees of freedom
117 * qinv = returned initial estimate
118 *
119 * vmax is cutoff above which degrees of freedom
120 * is treated as infinity.
121 */
122
qinv(double p,double c,double v)123 static double qinv(double p, double c, double v)
124 {
125 static const double p0 = 0.322232421088;
126 static const double q0 = 0.993484626060e-01;
127 static const double p1 = -1.0;
128 static const double q1 = 0.588581570495;
129 static const double p2 = -0.342242088547;
130 static const double q2 = 0.531103462366;
131 static const double p3 = -0.204231210125;
132 static const double q3 = 0.103537752850;
133 static const double p4 = -0.453642210148e-04;
134 static const double q4 = 0.38560700634e-02;
135 static const double c1 = 0.8832;
136 static const double c2 = 0.2368;
137 static const double c3 = 1.214;
138 static const double c4 = 1.208;
139 static const double c5 = 1.4142;
140 static const double vmax = 120.0;
141
142 double ps, q, t, yi;
143
144 ps = 0.5 - 0.5 * p;
145 yi = sqrt (log (1.0 / (ps * ps)));
146 t = yi + ((((yi * p4 + p3) * yi + p2) * yi + p1) * yi + p0)
147 / ((((yi * q4 + q3) * yi + q2) * yi + q1) * yi + q0);
148 if (v < vmax) t += (t * t * t + t) / v / 4.0;
149 q = c1 - c2 * t;
150 if (v < vmax) q += -c3 / v + c4 * t / v;
151 return t * (q * log (c - 1.0) + c5);
152 }
153
154 /*
155 * Copenhaver, Margaret Diponzio & Holland, Burt S.
156 * Multiple comparisons of simple effects in
157 * the two-way analysis of variance with fixed effects.
158 * Journal of Statistical Computation and Simulation,
159 * Vol.30, pp.1-15, 1988.
160 *
161 * Uses the secant method to find critical values.
162 *
163 * p = confidence level (1 - alpha)
164 * rr = no. of rows or groups
165 * cc = no. of columns or treatments
166 * df = degrees of freedom of error term
167 *
168 * ir(1) = error flag = 1 if wprob probability > 1
169 * ir(2) = error flag = 1 if ptukey probability > 1
170 * ir(3) = error flag = 1 if convergence not reached in 50 iterations
171 * = 2 if df < 2
172 *
173 * qtukey = returned critical value
174 *
175 * If the difference between successive iterates is less than eps,
176 * the search is terminated
177 */
178
179
qtukey(double p,double rr,double cc,double df,int lower_tail,int log_p)180 double qtukey(double p, double rr, double cc, double df,
181 int lower_tail, int log_p)
182 {
183 static const double eps = 0.0001;
184 const int maxiter = 50;
185
186 double ans = 0.0, valx0, valx1, x0, x1, xabs;
187 int iter;
188
189 if (isnan(p) || isnan(rr) || isnan(cc) || isnan(df)) {
190 /* ML_ERROR(ME_DOMAIN, "qtukey"); */
191 return p + rr + cc + df;
192 }
193
194 /* df must be > 1 ; there must be at least two values */
195 /* ^^
196 JMD: The comment says 1 but the code says 2.
197 Which is correct?
198 */
199 assert (df >= 2);
200 assert (rr >= 1);
201 assert (cc >= 2);
202
203
204 R_Q_P01_boundaries (p, 0, ML_POSINF);
205
206 p = R_DT_qIv(p); /* lower_tail,non-log "p" */
207
208 /* Initial value */
209
210 x0 = qinv(p, cc, df);
211
212 /* Find prob(value < x0) */
213
214 valx0 = ptukey(x0, rr, cc, df, /*LOWER*/TRUE, /*LOG_P*/FALSE) - p;
215
216 /* Find the second iterate and prob(value < x1). */
217 /* If the first iterate has probability value */
218 /* exceeding p then second iterate is 1 less than */
219 /* first iterate; otherwise it is 1 greater. */
220
221 if (valx0 > 0.0)
222 x1 = fmax2(0.0, x0 - 1.0);
223 else
224 x1 = x0 + 1.0;
225 valx1 = ptukey(x1, rr, cc, df, /*LOWER*/TRUE, /*LOG_P*/FALSE) - p;
226
227 /* Find new iterate */
228
229 for(iter=1 ; iter < maxiter ; iter++) {
230 ans = x1 - ((valx1 * (x1 - x0)) / (valx1 - valx0));
231 valx0 = valx1;
232
233 /* New iterate must be >= 0 */
234
235 x0 = x1;
236 if (ans < 0.0) {
237 ans = 0.0;
238 valx1 = -p;
239 }
240 /* Find prob(value < new iterate) */
241
242 valx1 = ptukey(ans, rr, cc, df, /*LOWER*/TRUE, /*LOG_P*/FALSE) - p;
243 x1 = ans;
244
245 /* If the difference between two successive */
246 /* iterates is less than eps, stop */
247
248 xabs = fabs(x1 - x0);
249 if (xabs < eps)
250 return ans;
251 }
252
253 /* The process did not converge in 'maxiter' iterations */
254 assert (0);
255 return ans;
256 }
257