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 #ifdef IEEE_754
  85         if (ISNAN(x) || ISNAN(y))
  86                 return x + y;
  87 #endif
  88         return (x < y) ? y : x;
  89 }
  90
  91
  92 #define R_Q_P01_boundaries(p, _LEFT_, _RIGHT_)          \
  93     if (log_p) {                                        \
  94       assert (p <= 0);                                  \
  95         if(p == 0) /* upper bound*/                     \
  96             return lower_tail ? _RIGHT_ : _LEFT_;       \
  97         if(p == ML_NEGINF)                              \
  98             return lower_tail ? _LEFT_ : _RIGHT_;       \
  99     }                                                   \
 100     else { /* !log_p */                                 \
 101       assert (p >= 0 && p <= 1);                        \
 102         if(p == 0)                                      \
 103             return lower_tail ? _LEFT_ : _RIGHT_;       \
 104         if(p == 1)                                      \
 105             return lower_tail ? _RIGHT_ : _LEFT_;       \
 106     }
 107
 108
 109 /* qinv() :
 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.
 115  *
 116  *        p = percentage point
 117  *        c = no. of columns or treatments
 118  *        v = degrees of freedom
 119  *        qinv = returned initial estimate
 120  *
 121  *      vmax is cutoff above which degrees of freedom
 122  *      is treated as infinity.
 123  */
 124
 125 static double qinv(double p, double c, double v)
 126 {
 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;
 143
 144     double ps, q, t, yi;
 145
 146     ps = 0.5 - 0.5 * p;
 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;
 151     q = c1 - c2 * t;
 152     if (v < vmax) q += -c3 / v + c4 * t / v;
 153     return t * (q * log (c - 1.0) + c5);
 154 }
 155
 156 /*
 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.
 162  *
 163  *  Uses the secant method to find critical values.
 164  *
 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
 169  *
 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
 173  *                     = 2 if df < 2
 174  *
 175  *  qtukey = returned critical value
 176  *
 177  *  If the difference between successive iterates is less than eps,
 178  *  the search is terminated
 179  */
 180
 181
 182 double qtukey(double p, double rr, double cc, double df,
 183               int lower_tail, int log_p)
 184 {
 185     static const double eps = 0.0001;
 186     const int maxiter = 50;
 187
 188     double ans = 0.0, valx0, valx1, x0, x1, xabs;
 189     int iter;
 190
 191 #ifdef IEEE_754
 192     if (ISNAN(p) || ISNAN(rr) || ISNAN(cc) || ISNAN(df)) {
 193         ML_ERROR(ME_DOMAIN, "qtukey");
 194         return p + rr + cc + df;
 195     }
 196 #endif
 197
 198     /* df must be > 1 ; there must be at least two values */
 199     assert (! (df < 2 || rr < 1 || cc < 2) );
 200
 201     R_Q_P01_boundaries (p, 0, ML_POSINF);
 202
 203     p = R_DT_qIv(p); /* lower_tail,non-log "p" */
 204
 205     /* Initial value */
 206
 207     x0 = qinv(p, cc, df);
 208
 209     /* Find prob(value < x0) */
 210
 211     valx0 = ptukey(x0, rr, cc, df, /*LOWER*/TRUE, /*LOG_P*/FALSE) - p;
 212
 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. */
 217
 218     if (valx0 > 0.0)
 219         x1 = fmax2(0.0, x0 - 1.0);
 220     else
 221         x1 = x0 + 1.0;
 222     valx1 = ptukey(x1, rr, cc, df, /*LOWER*/TRUE, /*LOG_P*/FALSE) - p;
 223
 224     /* Find new iterate */
 225
 226     for(iter=1 ; iter < maxiter ; iter++) {
 227         ans = x1 - ((valx1 * (x1 - x0)) / (valx1 - valx0));
 228         valx0 = valx1;
 229
 230         /* New iterate must be >= 0 */
 231
 232         x0 = x1;
 233         if (ans < 0.0) {
 234             ans = 0.0;
 235             valx1 = -p;
 236         }
 237         /* Find prob(value < new iterate) */
 238
 239         valx1 = ptukey(ans, rr, cc, df, /*LOWER*/TRUE, /*LOG_P*/FALSE) - p;
 240         x1 = ans;
 241
 242         /* If the difference between two successive */
 243         /* iterates is less than eps, stop */
 244
 245         xabs = fabs(x1 - x0);
 246         if (xabs < eps)
 247             return ans;
 248     }
 249
 250     /* The process did not converge in 'maxiter' iterations */
 251     assert (0);
 252     return ans;
 253 }