GRASS GIS 8 Programmer's Manual  8.4.0dev(2024)-df949ca623
as177.c
Go to the documentation of this file.
1 /*-Algorithm AS 177
2  * Expected Normal Order Statistics (Exact and Approximate),
3  * by J.P. Royston, 1982.
4  * Applied Statistics, 31(2):161-165.
5  *
6  * Translation to C by James Darrell McCauley, mccauley@ecn.purdue.edu.
7  *
8  * The functions Cdhc_nscor1() and Cdhc_nscor2() calculate the expected values
9  * of normal order statistics in exact or approximate form, respectively.
10  *
11  */
12 
13 #define NSTEP 721
14 #define H 0.025
15 
16 #include <math.h>
17 #include <stdio.h>
18 #include "local_proto.h"
19 
20 /* Local function prototypes */
21 static double Cdhc_alnfac(int j);
22 static double Cdhc_correc(int i, int n);
23 
24 /* exact calculation of normal scores */
25 void Cdhc_nscor1(double s[], int n, int n2, double work[], int *ifault)
26 {
27  double ani, c, c1, d, scor;
28  int i, j;
29 
30  *ifault = 3;
31  if (n2 != n / 2)
32  return;
33 
34  *ifault = 1;
35  if (n <= 1)
36  return;
37 
38  *ifault = 0;
39  if (n > 2000)
40  *ifault = 2;
41 
42  /* calculate the natural log of factorial(n) */
43  c1 = Cdhc_alnfac(n);
44  d = c1 - log((double)n);
45 
46  /* accumulate ordinates for calculation of integral for rankits */
47  for (i = 0; i < n2; ++i) {
48  ani = (double)n - i - 1;
49  c = c1 - d;
50  for (scor = 0.0, j = 0; j < NSTEP; ++j)
51  scor += work[0 * NSTEP + j] *
52  exp(work[1 * NSTEP + j] + work[2 * NSTEP + j] * i +
53  work[3 * NSTEP + j] * ani + c);
54  s[i] = scor * H;
55  d += log((double)(i + 1.0) / ani);
56  }
57 
58  return;
59 }
60 
61 void init(double work[])
62 {
63  double xstart = -9.0, pi2 = -0.918938533, xx;
64  int i;
65 
66  xx = xstart;
67 
68  /* set up arrays for calculation of integral */
69  for (i = 0; i < NSTEP; ++i) {
70  work[0 * NSTEP + i] = xx;
71  work[1 * NSTEP + i] = pi2 - xx * xx * 0.5;
72  work[2 * NSTEP + i] = log(Cdhc_alnorm(xx, 1));
73  work[3 * NSTEP + i] = log(Cdhc_alnorm(xx, 0));
74  xx = xstart + H * (i + 1.0);
75  }
76 
77  return;
78 }
79 
80 /*-Algorithm AS 177.2 Appl. Statist. (1982) Vol.31, No.2
81  * Natural logarithm of factorial for non-negative argument
82  */
83 static double Cdhc_alnfac(int j)
84 {
85  static double r[7] = {0.0, 0.0, 0.69314718056,
86  1.79175946923, 3.17805383035, 4.78749174278,
87  6.57925121101};
88  double w, z;
89 
90  if (j == 1)
91  return (double)1.0;
92  else if (j <= 7)
93  return r[j];
94 
95  w = (double)j + 1;
96  z = 1.0 / (w * w);
97 
98  return (w - 0.5) * log(w) - w + 0.918938522305 +
99  (((4.0 - 3.0 * z) * z - 14.0) * z + 420.0) / (5040.0 * w);
100 }
101 
102 /*-Algorithm AS 177.3 Appl. Statist. (1982) Vol.31, No.2
103  * Approximation for Rankits
104  */
105 void Cdhc_nscor2(double s[], int n, int n2, int *ifault)
106 {
107  static double eps[4] = {0.419885, 0.450536, 0.456936, 0.468488};
108  static double dl1[4] = {0.112063, 0.121770, 0.239299, 0.215159};
109  static double dl2[4] = {0.080122, 0.111348, -0.211867, -0.115049};
110  static double gam[4] = {0.474798, 0.469051, 0.208597, 0.259784};
111  static double lam[4] = {0.282765, 0.304856, 0.407708, 0.414093};
112  static double bb = -0.283833, d = -0.106136, b1 = 0.5641896;
113  double e1, e2, l1;
114  int i, k;
115 
116  *ifault = 3;
117  if (n2 != n / 2)
118  return;
119 
120  *ifault = 1;
121  if (n <= 1)
122  return;
123 
124  *ifault = 0;
125  if (n > 2000)
126  *ifault = 2;
127 
128  s[0] = b1;
129  if (n == 2)
130  return;
131 
132  /* calculate normal areas for 3 largest rankits */
133  k = (n2 < 3) ? n2 : 3;
134  for (i = 0; i < k; ++i) {
135  e1 = (1.0 + i - eps[i]) / (n + gam[i]);
136  e2 = pow(e1, lam[i]);
137  s[i] = e1 + e2 * (dl1[i] + e2 * dl2[i]) / n - Cdhc_correc(1 + i, n);
138  }
139 
140  if (n2 != k) {
141  /* calculate normal areas for remaining rankits */
142  for (i = 3; i < n2; ++i) {
143  l1 = lam[3] + bb / (1.0 + i + d);
144  e1 = (1.0 + i - eps[3]) / (n + gam[3]);
145  e2 = pow(e1, l1);
146  s[i] = e1 + e2 * (dl1[3] + e2 * dl2[3]) / n - Cdhc_correc(1 + i, n);
147  }
148  }
149 
150  /* convert normal tail areas to normal deviates */
151  for (i = 0; i < n2; ++i)
152  s[i] = -ppnd16(s[i]);
153 
154  return;
155 }
156 
157 /*-Algorithm AS 177.4 Appl. Statist. (1982) Vol.31, No.2
158  * Calculates Cdhc_correction for tail area of normal distribution
159  * corresponding to ith largest rankit in sample size n.
160  */
161 static double Cdhc_correc(int i, int n)
162 {
163  static double c1[7] = {9.5, 28.7, 1.9, 0.0, -7.0, -6.2, -1.6};
164  static double c2[7] = {-6.195e3, -9.569e3, -6.728e3, -17.614e3,
165  -8.278e3, -3.570e3, 1.075e3};
166  static double c3[7] = {9.338e4, 1.7516e5, 4.1040e5, 2.157e6,
167  2.376e6, 2.065e6, 2.065e6};
168  static double mic = 1.0e-6, c14 = 1.9e-5;
169  double an;
170 
171  if (i * n == 4)
172  return c14;
173 
174  if (i < 1 || i > 7)
175  return 0.0;
176  else if (i != 4 && n > 20)
177  return 0.0;
178  else if (i == 4 && n > 40)
179  return 0.0;
180 
181  /* else */
182  an = 1.0 / (double)(n * n);
183  return (c1[i - 1] + an * (c2[i - 1] + an * c3[i - 1])) * mic;
184 }
#define NSTEP
Definition: as177.c:13
void init(double work[])
Definition: as177.c:61
void Cdhc_nscor1(double s[], int n, int n2, double work[], int *ifault)
Definition: as177.c:25
void Cdhc_nscor2(double s[], int n, int n2, int *ifault)
Definition: as177.c:105
#define H
Definition: as177.c:14
double ppnd16(double p)
Definition: as241.c:86
double Cdhc_alnorm(double x, int upper)
Definition: as66.c:32
double r
Definition: r_raster.c:39