zip-reader: Read the whole central directory at .zip open time.
[pspp] / lib / tukey / qtukey.c
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
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
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
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 }