GRASS Programmer's Manual  6.4.2(2012)
as177.c
Go to the documentation of this file.
00001 
00002 /*-Algorithm AS 177
00003  * Expected Normal Order Statistics (Exact and Approximate),
00004  * by J.P. Royston, 1982.
00005  * Applied Statistics, 31(2):161-165.
00006  *
00007  * Translation to C by James Darrell McCauley, mccauley@ecn.purdue.edu.
00008  *
00009  * The functions nscor1() and nscor2() calculate the expected values of
00010  * normal order statistics in exact or approximate form, respectively.
00011  *
00012  */
00013 
00014 #define NSTEP 721
00015 #define H 0.025
00016 
00017 #include <math.h>
00018 #include <stdio.h>
00019 #include "local_proto.h"
00020 
00021 
00022 /* Local function prototypes */
00023 static double alnfac(int j);
00024 static double correc(int i, int n);
00025 
00026 
00027 /* exact calculation of normal scores */
00028 void nscor1(double s[], int n, int n2, double work[], int *ifault)
00029 {
00030     double ani, c, c1, d, scor;
00031     int i, j;
00032 
00033     *ifault = 3;
00034     if (n2 != n / 2)
00035         return;
00036 
00037     *ifault = 1;
00038     if (n <= 1)
00039         return;
00040 
00041     *ifault = 0;
00042     if (n > 2000)
00043         *ifault = 2;
00044 
00045     /* calculate the natural log of factorial(n) */
00046     c1 = alnfac(n);
00047     d = c1 - log((double)n);
00048 
00049     /* accumulate ordinates for calculation of integral for rankits */
00050     for (i = 0; i < n2; ++i) {
00051         ani = (double)n - i - 1;
00052         c = c1 - d;
00053         for (scor = 0.0, j = 0; j < NSTEP; ++j)
00054             scor += work[0 * NSTEP + j] *
00055                 exp(work[1 * NSTEP + j] + work[2 * NSTEP + j] * i
00056                     + work[3 * NSTEP + j] * ani + c);
00057         s[i] = scor * H;
00058         d += log((double)(i + 1.0) / ani);
00059     }
00060 
00061     return;
00062 }
00063 
00064 
00065 void init(double work[])
00066 {
00067     double xstart = -9.0, pi2 = -0.918938533, xx;
00068     int i;
00069 
00070     xx = xstart;
00071 
00072     /* set up arrays for calculation of integral */
00073     for (i = 0; i < NSTEP; ++i) {
00074         work[0 * NSTEP + i] = xx;
00075         work[1 * NSTEP + i] = pi2 - xx * xx * 0.5;
00076         work[2 * NSTEP + i] = log(alnorm(xx, 1));
00077         work[3 * NSTEP + i] = log(alnorm(xx, 0));
00078         xx = xstart + H * (i + 1.0);
00079     }
00080 
00081     return;
00082 }
00083 
00084 
00085 /*-Algorithm AS 177.2 Appl. Statist. (1982) Vol.31, No.2
00086  * Natural logarithm of factorial for non-negative argument
00087  */
00088 static double alnfac(int j)
00089 {
00090     static double r[7] = { 0.0, 0.0, 0.69314718056, 1.79175946923,
00091         3.17805383035, 4.78749174278, 6.57925121101
00092     };
00093     double w, z;
00094 
00095     if (j == 1)
00096         return (double)1.0;
00097     else if (j <= 7)
00098         return r[j];
00099 
00100     w = (double)j + 1;
00101     z = 1.0 / (w * w);
00102 
00103     return (w - 0.5) * log(w) - w + 0.918938522305 +
00104         (((4.0 - 3.0 * z) * z - 14.0) * z + 420.0) / (5040.0 * w);
00105 }
00106 
00107 
00108 /*-Algorithm AS 177.3 Appl. Statist. (1982) Vol.31, No.2
00109  * Approximation for Rankits
00110  */
00111 void nscor2(double s[], int n, int n2, int *ifault)
00112 {
00113     static double eps[4] = { 0.419885, 0.450536, 0.456936, 0.468488 };
00114     static double dl1[4] = { 0.112063, 0.121770, 0.239299, 0.215159 };
00115     static double dl2[4] = { 0.080122, 0.111348, -0.211867, -0.115049 };
00116     static double gam[4] = { 0.474798, 0.469051, 0.208597, 0.259784 };
00117     static double lam[4] = { 0.282765, 0.304856, 0.407708, 0.414093 };
00118     static double bb = -0.283833, d = -0.106136, b1 = 0.5641896;
00119     double e1, e2, l1;
00120     int i, k;
00121 
00122     *ifault = 3;
00123     if (n2 != n / 2)
00124         return;
00125 
00126     *ifault = 1;
00127     if (n <= 1)
00128         return;
00129 
00130     *ifault = 0;
00131     if (n > 2000)
00132         *ifault = 2;
00133 
00134     s[0] = b1;
00135     if (n == 2)
00136         return;
00137 
00138     /* calculate normal areas for 3 largest rankits */
00139     k = (n2 < 3) ? n2 : 3;
00140     for (i = 0; i < k; ++i) {
00141         e1 = (1.0 + i - eps[i]) / (n + gam[i]);
00142         e2 = pow(e1, lam[i]);
00143         s[i] = e1 + e2 * (dl1[i] + e2 * dl2[i]) / n - correc(1 + i, n);
00144     }
00145 
00146     if (n2 != k) {
00147         /* calculate normal areas for remaining rankits */
00148         for (i = 3; i < n2; ++i) {
00149             l1 = lam[3] + bb / (1.0 + i + d);
00150             e1 = (1.0 + i - eps[3]) / (n + gam[3]);
00151             e2 = pow(e1, l1);
00152             s[i] = e1 + e2 * (dl1[3] + e2 * dl2[3]) / n - correc(1 + i, n);
00153         }
00154     }
00155 
00156     /* convert normal tail areas to normal deviates */
00157     for (i = 0; i < n2; ++i)
00158         s[i] = -ppnd16(s[i]);
00159 
00160     return;
00161 }
00162 
00163 
00164 /*-Algorithm AS 177.4 Appl. Statist. (1982) Vol.31, No.2
00165  * Calculates correction for tail area of noraml distribution
00166  * corresponding to ith largest rankit in sample size n.
00167  */
00168 static double correc(int i, int n)
00169 {
00170     static double c1[7] = { 9.5, 28.7, 1.9, 0.0, -7.0, -6.2, -1.6 };
00171     static double c2[7] = { -6.195e3, -9.569e3, -6.728e3, -17.614e3,
00172         -8.278e3, -3.570e3, 1.075e3
00173     };
00174     static double c3[7] = { 9.338e4, 1.7516e5, 4.1040e5, 2.157e6,
00175         2.376e6, 2.065e6, 2.065e6
00176     };
00177     static double mic = 1.0e-6, c14 = 1.9e-5;
00178     double an;
00179 
00180     if (i * n == 4)
00181         return c14;
00182 
00183     if (i < 1 || i > 7)
00184         return 0.0;
00185     else if (i != 4 && n > 20)
00186         return 0.0;
00187     else if (i == 4 && n > 40)
00188         return 0.0;
00189 
00190     /* else */
00191     an = 1.0 / (double)(n * n);
00192     return (c1[i - 1] + an * (c2[i - 1] + an * c3[i - 1])) * mic;
00193 }
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines