GRASS Programmer's Manual
6.4.2(2012)
|
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 }