1 /*
2 * R : A Computer Language for Statistical Data Analysis
3 * Copyright (C) 2000-2016 The R Core Team
4 * Copyright (C) 2003 The R Foundation
5 * based on AS 89 (C) 1975 Royal Statistical Society
6 *
7 * This program is free software; you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation; either version 2 of the License, or
10 * (at your option) any later version.
11 *
12 * This program is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with this program; if not, a copy is available at
19 * https://www.R-project.org/Licenses/
20 *
21 */
22
23 #include <math.h>
24 #include <Rmath.h>
25
26 /* Was
27 double precision function prho(n, is, ifault)
28
29 Changed to subroutine by KH to allow for .Fortran interfacing.
30 Also, change `prho' argument in the code to `pv'.
31 And, fix a bug.
32
33 From R ver. 1.1.x [March, 2000] by MM:
34 - Translate Fortran to C
35 - use pnorm() instead of less precise alnorm().
36 - new argument lower_tail --> potentially increased precision in extreme cases.
37 */
38 void
prho(int n,double is,double * pv,int ifault,int lower_tail)39 prho(int n, double is, double *pv, int ifault, int lower_tail)
40 {
41 /* Algorithm AS 89 Appl. Statist. (1975) Vol.24, No. 3, P377.
42
43 To evaluate the probability Pr[ S >= is ]
44 {or Pr [ S < is] if(lower_tail) }, where
45
46 S = (n^3 - n) * (1-R)/6,
47 is = (n^3 - n) * (1-r)/6,
48 R,r = Spearman's rho (r.v. and observed), and n >= 2
49 */
50
51 /* Edgeworth coefficients : */
52 const double
53 c1 = .2274,
54 c2 = .2531,
55 c3 = .1745,
56 c4 = .0758,
57 c5 = .1033,
58 c6 = .3932,
59 c7 = .0879,
60 c8 = .0151,
61 c9 = .0072,
62 c10= .0831,
63 c11= .0131,
64 c12= 4.6e-4;
65
66 /* Local variables */
67 double b, u, x, y, n3;/*, js */
68
69 #define n_small 9
70 /* originally: n_small = 6 (speed!);
71 * even n_small = 10 (and n = 10) needs quite a bit longer than the approx!
72 * larger than 12 ==> integer overflow in nfac and (probably) ifr
73 */
74 int l[n_small];
75 int nfac, i, m, mt, ifr, ise, n1;
76
77 /* Test admissibility of arguments and initialize */
78 *pv = lower_tail ? 0. : 1.;
79 if (n <= 1) { ifault = 1; return; }
80
81 ifault = 0;
82 if (is <= 0.) return;/* with p = 1 */
83
84 n3 = (double)n;
85 n3 *= (n3 * n3 - 1.) / 3.;/* = (n^3 - n)/3 */
86 if (is > n3) { /* larger than maximal value */
87 *pv = 1 - *pv; return;
88 }
89 /* NOT rounding to even anymore: with ties, S, may even be non-integer!
90 * js = is;
91 * if(fmod(js, 2.) != 0.) ++js;
92 */
93 if (n <= n_small) { /* 2 <= n <= n_small :
94 * Exact evaluation of probability */
95 nfac = 1.;
96 for (i = 1; i <= n; ++i) {
97 nfac *= i;
98 l[i - 1] = i;
99 }
100 /* KH mod next line: was `!=' in the code but `.eq.' in the paper */
101 if (is == n3) {
102 ifr = 1;
103 }
104 else {
105 ifr = 0;
106 for (m = 0; m < nfac; ++m) {
107 ise = 0;
108 for (i = 0; i < n; ++i) {
109 n1 = i + 1 - l[i];
110 ise += n1 * n1;
111 }
112 if (is <= ise)
113 ++ifr;
114
115 n1 = n;
116 do {
117 mt = l[0];
118 for (i = 1; i < n1; ++i)
119 l[i - 1] = l[i];
120 --n1;
121 l[n1] = mt;
122 } while (mt == n1+1 && n1 > 1);
123 }
124 }
125 *pv = (lower_tail ? nfac-ifr : ifr) / (double) nfac;
126 } /* exact for n <= n_small */
127
128 else { /* n >= 7 : Evaluation by Edgeworth series expansion */
129
130 y = (double) n;
131 b = 1 / y;
132 x = (6. * (is - 1) * b / (y * y - 1) - 1) * sqrt(y - 1);
133 /* = rho * sqrt(n-1) == rho / sqrt(var(rho)) ~ (0,1) */
134 y = x * x;
135 u = x * b * (c1 + b * (c2 + c3 * b) +
136 y * (-c4 + b * (c5 + c6 * b) -
137 y * b * (c7 + c8 * b -
138 y * (c9 - c10 * b + y * b * (c11 - c12 * y))
139 )));
140 y = u / exp(y / 2.);
141 *pv = (lower_tail ? -y : y) +
142 pnorm(x, 0., 1., lower_tail, /*log_p = */FALSE);
143 /* above was call to alnorm() [algorithm AS 66] */
144 if (*pv < 0) *pv = 0.;
145 if (*pv > 1) *pv = 1.;
146 }
147 return;
148 } /* prho */
149
150 #include <Rinternals.h>
pRho(SEXP q,SEXP sn,SEXP lower)151 SEXP pRho(SEXP q, SEXP sn, SEXP lower)
152 {
153 double s = asReal(q), p;
154 int n = asInteger(sn), ltail = asInteger(lower), ifault = 0;
155 prho(n, s, &p, ifault, ltail);
156 return ScalarReal(p);
157 }
158