1 /* PSPP - a program for statistical analysis.
2 Copyright (C) 2011 Free Software Foundation, Inc.
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.
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.
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/>.
18 /* This file is taken from the R project source code, and modified.
19 The original copyright notice is reproduced below: */
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
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.
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.
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/
44 * double qtukey(p, rr, cc, df, lower_tail, log_p);
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.
52 * The algorithm is based on that of the reference.
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.
73 #define ML_POSINF (1.0 / 0.0)
74 #define ML_NEGINF (-1.0 / 0.0)
76 #define R_D_Lval(p) (lower_tail ? (p) : (0.5 - (p) + 0.5)) /* p */
78 #define R_DT_qIv(p) (log_p ? (lower_tail ? exp(p) : - expm1(p)) \
82 static double fmax2(double x, double y)
85 if (ISNAN(x) || ISNAN(y))
88 return (x < y) ? y : x;
92 #define R_Q_P01_boundaries(p, _LEFT_, _RIGHT_) \
95 if(p == 0) /* upper bound*/ \
96 return lower_tail ? _RIGHT_ : _LEFT_; \
98 return lower_tail ? _LEFT_ : _RIGHT_; \
100 else { /* !log_p */ \
101 assert (p >= 0 && p <= 1); \
103 return lower_tail ? _LEFT_ : _RIGHT_; \
105 return lower_tail ? _RIGHT_ : _LEFT_; \
110 * this function finds percentage point of the studentized range
111 * which is used as initial estimate for the secant method.
112 * function is adapted from portion of algorithm as 70
113 * from applied statistics (1974) ,vol. 23, no. 1
114 * by odeh, r. e. and evans, j. o.
116 * p = percentage point
117 * c = no. of columns or treatments
118 * v = degrees of freedom
119 * qinv = returned initial estimate
121 * vmax is cutoff above which degrees of freedom
122 * is treated as infinity.
125 static double qinv(double p, double c, double v)
127 static const double p0 = 0.322232421088;
128 static const double q0 = 0.993484626060e-01;
129 static const double p1 = -1.0;
130 static const double q1 = 0.588581570495;
131 static const double p2 = -0.342242088547;
132 static const double q2 = 0.531103462366;
133 static const double p3 = -0.204231210125;
134 static const double q3 = 0.103537752850;
135 static const double p4 = -0.453642210148e-04;
136 static const double q4 = 0.38560700634e-02;
137 static const double c1 = 0.8832;
138 static const double c2 = 0.2368;
139 static const double c3 = 1.214;
140 static const double c4 = 1.208;
141 static const double c5 = 1.4142;
142 static const double vmax = 120.0;
147 yi = sqrt (log (1.0 / (ps * ps)));
148 t = yi + (((( yi * p4 + p3) * yi + p2) * yi + p1) * yi + p0)
149 / (((( yi * q4 + q3) * yi + q2) * yi + q1) * yi + q0);
150 if (v < vmax) t += (t * t * t + t) / v / 4.0;
152 if (v < vmax) q += -c3 / v + c4 * t / v;
153 return t * (q * log (c - 1.0) + c5);
157 * Copenhaver, Margaret Diponzio & Holland, Burt S.
158 * Multiple comparisons of simple effects in
159 * the two-way analysis of variance with fixed effects.
160 * Journal of Statistical Computation and Simulation,
161 * Vol.30, pp.1-15, 1988.
163 * Uses the secant method to find critical values.
165 * p = confidence level (1 - alpha)
166 * rr = no. of rows or groups
167 * cc = no. of columns or treatments
168 * df = degrees of freedom of error term
170 * ir(1) = error flag = 1 if wprob probability > 1
171 * ir(2) = error flag = 1 if ptukey probability > 1
172 * ir(3) = error flag = 1 if convergence not reached in 50 iterations
175 * qtukey = returned critical value
177 * If the difference between successive iterates is less than eps,
178 * the search is terminated
182 double qtukey(double p, double rr, double cc, double df,
183 int lower_tail, int log_p)
185 static const double eps = 0.0001;
186 const int maxiter = 50;
188 double ans = 0.0, valx0, valx1, x0, x1, xabs;
192 if (ISNAN(p) || ISNAN(rr) || ISNAN(cc) || ISNAN(df)) {
193 ML_ERROR(ME_DOMAIN, "qtukey");
194 return p + rr + cc + df;
198 /* df must be > 1 ; there must be at least two values */
199 assert (! (df < 2 || rr < 1 || cc < 2) );
201 R_Q_P01_boundaries (p, 0, ML_POSINF);
203 p = R_DT_qIv(p); /* lower_tail,non-log "p" */
207 x0 = qinv(p, cc, df);
209 /* Find prob(value < x0) */
211 valx0 = ptukey(x0, rr, cc, df, /*LOWER*/TRUE, /*LOG_P*/FALSE) - p;
213 /* Find the second iterate and prob(value < x1). */
214 /* If the first iterate has probability value */
215 /* exceeding p then second iterate is 1 less than */
216 /* first iterate; otherwise it is 1 greater. */
219 x1 = fmax2(0.0, x0 - 1.0);
222 valx1 = ptukey(x1, rr, cc, df, /*LOWER*/TRUE, /*LOG_P*/FALSE) - p;
224 /* Find new iterate */
226 for(iter=1 ; iter < maxiter ; iter++) {
227 ans = x1 - ((valx1 * (x1 - x0)) / (valx1 - valx0));
230 /* New iterate must be >= 0 */
237 /* Find prob(value < new iterate) */
239 valx1 = ptukey(ans, rr, cc, df, /*LOWER*/TRUE, /*LOG_P*/FALSE) - p;
242 /* If the difference between two successive */
243 /* iterates is less than eps, stop */
245 xabs = fabs(x1 - x0);
250 /* The process did not converge in 'maxiter' iterations */