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