1 /*
2 * Algorithm AS 275 Appl.Statist. (1992), vol.41, no.2
3 * original (C) 1992 Royal Statistical Society
4 *
5 * Computes the noncentral chi-squared distribution function with
6 * positive real degrees of freedom df and nonnegative noncentrality
7 * parameter ncp. pnchisq_raw is based on
8 *
9 * Ding, C. G. (1992)
10 * Algorithm AS275: Computing the non-central chi-squared
11 * distribution function. Appl.Statist., 41, 478-482.
12
13 * Other parts
14 * Copyright (C) 2000-2012 The R Core Team
15 * Copyright (C) 2003-2009 The R Foundation
16 */
17
18 #include "nmath.h"
19 #include "dpq.h"
20
21 /*----------- DEBUGGING -------------
22 *
23 * make CFLAGS='-DDEBUG_pnch ....'
24
25 * -- Feb.6, 2000 (R pre0.99); M.Maechler: still have
26 * bad precision & non-convergence in some cases (x ~= f, both LARGE)
27 */
28
29 #ifdef HAVE_LONG_DOUBLE
30 # define EXP expl
31 # define FABS fabsl
32 #else
33 # define EXP exp
34 # define FABS fabs
35 #endif
36
pnchisq(double x,double df,double ncp,int lower_tail,int log_p)37 double pnchisq(double x, double df, double ncp, int lower_tail, int log_p)
38 {
39 double ans;
40 #ifdef IEEE_754
41 if (ISNAN(x) || ISNAN(df) || ISNAN(ncp))
42 return x + df + ncp;
43 if (!R_FINITE(df) || !R_FINITE(ncp))
44 ML_ERR_return_NAN;
45 #endif
46
47 if (df < 0. || ncp < 0.) ML_ERR_return_NAN;
48
49 ans = pnchisq_raw(x, df, ncp, 1e-12, 8*DBL_EPSILON, 1000000, lower_tail);
50 if(ncp >= 80) {
51 if(lower_tail) {
52 ans = fmin2(ans, 1.0); /* e.g., pchisq(555, 1.01, ncp = 80) */
53 } else { /* !lower_tail */
54 /* since we computed the other tail cancellation is likely */
55 if(ans < 1e-10) ML_ERROR(ME_PRECISION, "pnchisq");
56 ans = fmax2(ans, 0.0); /* Precaution PR#7099 */
57 }
58 }
59 if (!log_p) return ans;
60 /* if ans is near one, we can do better using the other tail */
61 if (ncp >= 80 || ans < 1 - 1e-8) return log(ans);
62 ans = pnchisq_raw(x, df, ncp, 1e-12, 8*DBL_EPSILON, 1000000, !lower_tail);
63 return log1p(-ans);
64 }
65
66 double attribute_hidden
pnchisq_raw(double x,double f,double theta,double errmax,double reltol,int itrmax,Rboolean lower_tail)67 pnchisq_raw(double x, double f, double theta,
68 double errmax, double reltol, int itrmax, Rboolean lower_tail)
69 {
70 double lam, x2, f2, term, bound, f_x_2n, f_2n;
71 double l_lam = -1., l_x = -1.; /* initialized for -Wall */
72 int n;
73 Rboolean lamSml, tSml, is_r, is_b, is_it;
74 LDOUBLE ans, u, v, t, lt, lu =-1;
75
76 static const double _dbl_min_exp = M_LN2 * DBL_MIN_EXP;
77 /*= -708.3964 for IEEE double precision */
78
79 if (x <= 0.) {
80 if(x == 0. && f == 0.)
81 return lower_tail ? exp(-0.5*theta) : -expm1(-0.5*theta);
82 /* x < 0 or {x==0, f > 0} */
83 return lower_tail ? 0. : 1.;
84 }
85 if(!R_FINITE(x)) return lower_tail ? 1. : 0.;
86
87 /* This is principally for use from qnchisq */
88 #ifndef MATHLIB_STANDALONE
89 R_CheckUserInterrupt();
90 #endif
91
92 if(theta < 80) { /* use 110 for Inf, as ppois(110, 80/2, lower.tail=FALSE) is 2e-20 */
93 LDOUBLE sum = 0, sum2 = 0, lambda = 0.5*theta,
94 pr = EXP(-lambda); // does this need a feature test?
95 double ans;
96 int i;
97 /* we need to renormalize here: the result could be very close to 1 */
98 for(i = 0; i < 110; pr *= lambda/++i) {
99 sum2 += pr;
100 sum += pr * pchisq(x, f+2*i, lower_tail, FALSE);
101 if (sum2 >= 1-1e-15) break;
102 }
103 ans = (double) (sum/sum2);
104 return ans;
105 }
106
107
108 #ifdef DEBUG_pnch
109 REprintf("pnchisq(x=%g, f=%g, theta=%g): ",x,f,theta);
110 #endif
111 lam = .5 * theta;
112 lamSml = (-lam < _dbl_min_exp);
113 if(lamSml) {
114 /* MATHLIB_ERROR(
115 "non centrality parameter (= %g) too large for current algorithm",
116 theta) */
117 u = 0;
118 lu = -lam;/* == ln(u) */
119 l_lam = log(lam);
120 } else {
121 u = exp(-lam);
122 }
123
124 /* evaluate the first term */
125 v = u;
126 x2 = .5 * x;
127 f2 = .5 * f;
128 f_x_2n = f - x;
129
130 #ifdef DEBUG_pnch
131 REprintf("-- v=exp(-th/2)=%g, x/2= %g, f/2= %g\n",v,x2,f2);
132 #endif
133
134 if(f2 * DBL_EPSILON > 0.125 && /* very large f and x ~= f: probably needs */
135 FABS(t = x2 - f2) < /* another algorithm anyway */
136 sqrt(DBL_EPSILON) * f2) {
137 /* evade cancellation error */
138 /* t = exp((1 - t)*(2 - t/(f2 + 1))) / sqrt(2*M_PI*(f2 + 1));*/
139 lt = (1 - t)*(2 - t/(f2 + 1)) - 0.5 * log(2*M_PI*(f2 + 1));
140 #ifdef DEBUG_pnch
141 REprintf(" (case I) ==> ");
142 #endif
143 }
144 else {
145 /* Usual case 2: careful not to overflow .. : */
146 lt = f2*log(x2) -x2 - lgammafn(f2 + 1);
147 }
148 #ifdef DEBUG_pnch
149 REprintf(" lt= %g", lt);
150 #endif
151
152 tSml = (lt < _dbl_min_exp);
153 if(tSml) {
154 if (x > f + theta + 5* sqrt( 2*(f + 2*theta))) {
155 /* x > E[X] + 5* sigma(X) */
156 return lower_tail ? 1. : 0.; /* FIXME: We could be more accurate than 0. */
157 } /* else */
158 l_x = log(x);
159 ans = term = 0.; t = 0;
160 }
161 else {
162 t = EXP(lt);
163 #ifdef DEBUG_pnch
164 REprintf(", t=exp(lt)= %g\n", t);
165 #endif
166 ans = term = (double) (v * t);
167 }
168
169 for (n = 1, f_2n = f + 2., f_x_2n += 2.; ; n++, f_2n += 2, f_x_2n += 2) {
170 #ifdef DEBUG_pnch
171 REprintf("\n _OL_: n=%d",n);
172 #endif
173 #ifndef MATHLIB_STANDALONE
174 if(n % 1000) R_CheckUserInterrupt();
175 #endif
176 /* f_2n === f + 2*n
177 * f_x_2n === f - x + 2*n > 0 <==> (f+2n) > x */
178 if (f_x_2n > 0) {
179 /* find the error bound and check for convergence */
180
181 bound = (double) (t * x / f_x_2n);
182 #ifdef DEBUG_pnch
183 REprintf("\n L10: n=%d; term= %g; bound= %g",n,term,bound);
184 #endif
185 is_r = is_it = FALSE;
186 /* convergence only if BOTH absolute and relative error < 'bnd' */
187 if (((is_b = (bound <= errmax)) &&
188 (is_r = (term <= reltol * ans))) || (is_it = (n > itrmax)))
189 {
190 #ifdef DEBUG_pnch
191 REprintf("BREAK n=%d %s; bound= %g %s, rel.err= %g %s\n",
192 n, (is_it ? "> itrmax" : ""),
193 bound, (is_b ? "<= errmax" : ""),
194 term/ans, (is_r ? "<= reltol" : ""));
195 #endif
196 break; /* out completely */
197 }
198
199 }
200
201 /* evaluate the next term of the */
202 /* expansion and then the partial sum */
203
204 if(lamSml) {
205 lu += l_lam - log(n); /* u = u* lam / n */
206 if(lu >= _dbl_min_exp) {
207 /* no underflow anymore ==> change regime */
208 #ifdef DEBUG_pnch
209 REprintf(" n=%d; nomore underflow in u = exp(lu) ==> change\n",
210 n);
211 #endif
212 v = u = EXP(lu); /* the first non-0 'u' */
213 lamSml = FALSE;
214 }
215 } else {
216 u *= lam / n;
217 v += u;
218 }
219 if(tSml) {
220 lt += l_x - log(f_2n);/* t <- t * (x / f2n) */
221 if(lt >= _dbl_min_exp) {
222 /* no underflow anymore ==> change regime */
223 #ifdef DEBUG_pnch
224 REprintf(" n=%d; nomore underflow in t = exp(lt) ==> change\n",
225 n);
226 #endif
227 t = EXP(lt); /* the first non-0 't' */
228 tSml = FALSE;
229 }
230 } else {
231 t *= x / f_2n;
232 }
233 if(!lamSml && !tSml) {
234 term = (double) (v * t);
235 ans += term;
236 }
237
238 } /* for(n ...) */
239
240 if (is_it) {
241 MATHLIB_WARNING2(_("pnchisq(x=%g, ..): not converged in %d iter."),
242 x, itrmax);
243 }
244 #ifdef DEBUG_pnch
245 REprintf("\n == L_End: n=%d; term= %g; bound=%g\n",n,term,bound);
246 #endif
247 return (double) (lower_tail ? ans : 1 - ans);
248 }
249