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 }