1 /*
2 * R : A Computer Language for Statistical Data Analysis
3 * Copyright (C) 2000-1 The R Core Team.
4 *
5 * Based on Applied Statistics algorithm AS177
6 * (C) Royal Statistical Society 1982
7 *
8 * This program is free software; you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation; either version 2 of the License, or
11 * (at your option) any later version.
12 *
13 * This program is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public License
19 * along with this program; if not, a copy is available at
20 * https://www.R-project.org/Licenses/
21 */
22
23
24 /* nscor.f -- translated by f2c (version 19980913).
25 * ------- and produced by f2c-clean,v 1.8 --- and hand polished: M.Maechler
26 */
27
28 #include <Rmath.h>
29
30 #define nstep 721 /* = nrow(work) */
31
nscor1(double * s,int * n,int * n2,double * work,int * ifault)32 void nscor1(double *s, int *n, int *n2,
33 double *work, int *ifault)
34 {
35 /* Algorithm AS 177 Appl. Statist. (1982) Vol. 31, No. 2, 161-165
36
37 Exact calculation of Normal Scores
38 */
39
40 /* Initialized data */
41
42 const double hh = .025;
43
44 /* Local variables */
45 double scor, c, an, ai1, ani;
46 int i, j, i1, ni;
47
48 /* Parameter adjustments */
49 --s;
50 work -= 5;
51
52 /* Function Body */
53
54 if (*n2 != *n / 2) {
55 *ifault = 3; return;
56 }
57 if (*n <= 1) {
58 *ifault = 1; return;
59 }
60 *ifault = 0;
61 if (*n > 2000) {
62 *ifault = 2;/* non-fatal potential accuracy loss */
63 }
64
65 an = (double) (*n);
66 c = log(an);
67
68 /* Accumulate ordinates for calculation of integral for rankits */
69
70 for (i = 1; i <= *n2; ++i) {
71 i1 = i - 1;
72 ni = *n - i;
73 ai1 = (double) i1;
74 ani = (double) ni;
75 scor = 0.;
76 for (j = 1; j <= nstep; ++j) {
77 scor += exp( work[(j << 2) + 2] + ai1 * work[(j << 2) + 3] +
78 ani * work[(j << 2) + 4] + c) * work[(j << 2) + 1];
79 }
80 s[i] = scor * hh;
81 c += log(ani / (double) i);
82 }
83 return;
84 } /* nscor1 */
85
86
87
init(double * work)88 void init(double *work)
89 {
90 /* Initialize the work[] table for nscor1(.).
91
92 Algorithm AS 177.1 Appl. Statist. (1982) Vol. 31, No. 2
93
94 Now calling R's pnorm() instead of
95 external alnorm { = Appl.Stat. AS 66 }
96 */
97
98 const double xstart = -9.;
99 const double hh = .025;
100 const double pi2 = -.918938533;
101 const double half = .5;
102
103 int i;
104 double xx;
105
106 /* Parameter adjustments */
107 work -= 5;
108
109 xx = xstart;
110
111 /* Set up arrays for calculation of integral */
112
113 for (i = 1; i <= nstep; ++i) {
114 work[(i << 2) + 1] = xx;
115 work[(i << 2) + 2] = pi2 - xx * xx * half;
116 /* upper & lower tail */
117 /* had log(alnorm_(&xx, UPPER)) & log(alnorm_(&xx, LOWER)) :*/
118 work[(i << 2) + 3] = pnorm(xx,0.,1., 0, /*log_p = */1);
119 work[(i << 2) + 4] = pnorm(xx,0.,1., 1, /*log_p = */1);
120 xx = xstart + (double) i * hh;
121 }
122 return;
123 } /* init */
124
125
126 static double correc(int, int);
127
nscor2(float * s,int * n,int * n2,int * ier)128 void nscor2(float *s, int *n, int *n2, int *ier)
129 {
130
131 /* algorithm as 177.3, applied statistics, v.31, 161-165, 1982.
132
133 calculates approximate expected values of normal order statistics.
134 claimed accuracy is 0.0001, though usually accurate to 5-6 dec.
135
136 *** N.B. This routine is NOT in double precision ***
137
138 Arguments:
139
140 s(n2) = output, the first n2 expected values.
141 n = input, the sample size.
142 n2 = input, the number of order statistics required; must
143 be <= n/2.
144 ier = output, error indicator
145 = 0 if no error detected
146 = 1 if n <= 1.
147 = 2 if n > 2000, in which case the order statistics
148 are still calculated, but may be inaccurate.
149 = 3 if n2 > n/2 (n.b. this differs from the
150 published algorithm which returns an error
151 if n2 is not equal to n/2.)
152
153 Calls qnorm() [from R] which is an improvement of
154 ppnd = applied statistics algorithm 111.
155 An alternative is ppnd7 in algorithm AS 241.
156 */
157
158 /* Initialized data */
159
160 const float
161 eps[4] = { .419885f,.450536f, .456936f, .468488f },
162 dl1[4] = { .112063f,.12177f, .239299f, .215159f },
163 dl2[4] = { .080122f,.111348f,-.211867f,-.115049f },
164 gam[4] = { .474798f,.469051f, .208597f, .259784f },
165 lam[4] = { .282765f,.304856f,.407708f,.414093f };
166 const float bb = -.283833f;
167 const float d = -.106136f;
168 const float b1 = .5641896f;
169
170
171 /* Local variables */
172 int i, k;
173 float e1, e2, ai, an;
174
175 /* input parameter checks. */
176
177 if (*n2 > *n / 2) {
178 *ier = 3; return;
179 }
180 if (*n <= 1) {
181 *ier = 1; return;
182 }
183 *ier = 0;
184 if (*n > 2000) {
185 *ier = 2;
186 }
187 s[0] = b1;
188 if (*n == 2) {
189 return;
190 }
191
192 /* calculate normal tail areas for first 3 order statistics. */
193
194 an = (float) (*n);
195 k = 3;
196 if (*n2 < k)
197 k = *n2;
198 /* k := min(3, *n2) */
199 for (i = 0; i < k; ++i) {
200 ai = (float) i+1;
201 e1 = (ai - eps[i]) / (an + gam[i]);
202 e2 = pow((double) e1, (double) lam[i]);
203 s[i] = e1 + e2 * (dl1[i] + e2 * dl2[i]) / an - correc(i+1, *n);
204 }
205 if (*n2 > k) {
206
207 /* calculate normal areas for other cases. */
208
209 for (i = 4-1; i < *n2; ++i) {
210 ai = (float) i+1;
211 e1 = (ai - eps[3]) / (an + gam[3]);
212 e2 = pow((double) e1, (double) lam[3] + bb / (ai + d));
213 s[i] = e1 + e2 * (dl1[3] + e2 * dl2[3]) / an - correc(i+1, *n);
214 }
215 }
216 /* convert tail areas to normal deviates. */
217
218 for (i = 0; i < *n2; ++i)
219 s[i] = - qnorm(s[i], 0., 1., 1, 0);
220
221 return;
222 } /* nscor2 */
223
224
correc(int i,int n)225 static double correc(int i, int n)
226 {
227 /* calculates correction for tail area of the i-th largest of n
228 order statistics. */
229
230 const float
231 c1[7] = { 9.5f,28.7f,1.9f,0.f,-7.f,-6.2f,-1.6f },
232 c2[7] = { -6195.f,-9569.f,-6728.f,-17614.f,-8278.f,-3570.f, 1075.f },
233 c3[7] = { 93380.f,175160.f,410400.f,2157600.f,2.376e6f,
234 2.065e6f,2.065e6f };
235 const float mic = 1e-6f;
236 const float c14 = 1.9e-5f;
237
238 double an;
239
240 if (i * n == 4) return c14;
241 if (i < 1 || i > 7) return 0;
242 if (i != 4 && n > 20) return 0;
243 if (i == 4 && n > 40) return 0;
244 /* else : */
245 an = (double) n;
246 an = 1. / (an * an);
247 i--;
248 return((c1[i] + an * (c2[i] + an * c3[i])) * mic);
249 } /* correc */
250