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