1 /*
2 * Mathlib : A C Library of Special Functions
3 * Copyright (C) 1998-2015 The R Core Team
4 * based on AS243 (C) 1989 Royal Statistical Society
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 2 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program; if not, a copy is available at
18 * https://www.R-project.org/Licenses/
19 */
20
21 /* Algorithm AS 243 Lenth,R.V. (1989). Appl. Statist., Vol.38, 185-189.
22 * ----------------
23 * Cumulative probability at t of the non-central t-distribution
24 * with df degrees of freedom (may be fractional) and non-centrality
25 * parameter delta.
26 *
27 * NOTE
28 *
29 * Requires the following auxiliary routines:
30 *
31 * lgammafn(x) - log gamma function
32 * pbeta(x, a, b) - incomplete beta function
33 * pnorm(x) - normal distribution function
34 *
35 * CONSTANTS
36 *
37 * M_SQRT_2dPI = 1/ {gamma(1.5) * sqrt(2)} = sqrt(2 / pi)
38 * M_LN_SQRT_PI = ln(sqrt(pi)) = ln(pi)/2
39 */
40
41 #include "nmath.h"
42 #include "dpq.h"
43
44 /*----------- DEBUGGING -------------
45 *
46 * make CFLAGS='-DDEBUG_pnt -g'
47
48 * -- Feb.3, 1999; M.Maechler:
49 - For 't > ncp > 20' (or so) the result is completely WRONG! <== no longer true
50 - but for ncp > 100
51 */
52
pnt(double t,double df,double ncp,int lower_tail,int log_p)53 double pnt(double t, double df, double ncp, int lower_tail, int log_p)
54 {
55 double albeta, a, b, del, errbd, lambda, rxb, tt, x;
56 LDOUBLE geven, godd, p, q, s, tnc, xeven, xodd;
57 int it, negdel;
58
59 /* note - itrmax and errmax may be changed to suit one's needs. */
60
61 const int itrmax = 1000;
62 const static double errmax = 1.e-12;
63
64 if (df <= 0.0) ML_WARN_return_NAN;
65 if(ncp == 0.0) return pt(t, df, lower_tail, log_p);
66
67 if(!R_FINITE(t))
68 return (t < 0) ? R_DT_0 : R_DT_1;
69 if (t >= 0.) {
70 negdel = FALSE; tt = t; del = ncp;
71 }
72 else {
73 /* We deal quickly with left tail if extreme,
74 since pt(q, df, ncp) <= pt(0, df, ncp) = \Phi(-ncp) */
75 if (ncp > 40 && (!log_p || !lower_tail)) return R_DT_0;
76 negdel = TRUE; tt = -t; del = -ncp;
77 }
78
79 if (df > 4e5 || del*del > 2*M_LN2*(-(DBL_MIN_EXP))) {
80 /*-- 2nd part: if del > 37.62, then p=0 below
81 FIXME: test should depend on `df', `tt' AND `del' ! */
82 /* Approx. from Abramowitz & Stegun 26.7.10 (p.949) */
83 s = 1./(4.*df);
84
85 return pnorm((double)(tt*(1. - s)), del,
86 sqrt((double) (1. + tt*tt*2.*s)),
87 lower_tail != negdel, log_p);
88 }
89
90 /* initialize twin series */
91 /* Guenther, J. (1978). Statist. Computn. Simuln. vol.6, 199. */
92
93 x = t * t;
94 rxb = df/(x + df);/* := (1 - x) {x below} -- but more accurately */
95 x = x / (x + df);/* in [0,1) */
96 #ifdef DEBUG_pnt
97 REprintf("pnt(t=%7g, df=%7g, ncp=%7g) ==> x= %10g:",t,df,ncp, x);
98 #endif
99 if (x > 0.) {/* <==> t != 0 */
100 lambda = del * del;
101 p = .5 * exp(-.5 * lambda);
102 #ifdef DEBUG_pnt
103 REprintf("\t p=%10Lg\n",p);
104 #endif
105 if(p == 0.) { /* underflow! */
106
107 /*========== really use an other algorithm for this case !!! */
108 ML_WARNING(ME_UNDERFLOW, "pnt");
109 ML_WARNING(ME_RANGE, "pnt"); /* |ncp| too large */
110 return R_DT_0;
111 }
112 #ifdef DEBUG_pnt
113 REprintf("it 1e5*(godd, geven)| p q s"
114 /* 1.3 1..4..7.9 1..4..7.9|1..4..7.901 1..4..7.901 1..4..7.901 */
115 " pnt(*) errbd\n");
116 /* 1..4..7..0..34 1..4..7.9*/
117 #endif
118 q = M_SQRT_2dPI * p * del;
119 s = .5 - p;
120 /* s = 0.5 - p = 0.5*(1 - exp(-.5 L)) = -0.5*expm1(-.5 L)) */
121 if(s < 1e-7)
122 s = -0.5 * expm1(-0.5 * lambda);
123 a = .5;
124 b = .5 * df;
125 /* rxb = (1 - x) ^ b [ ~= 1 - b*x for tiny x --> see 'xeven' below]
126 * where '(1 - x)' =: rxb {accurately!} above */
127 rxb = pow(rxb, b);
128 albeta = M_LN_SQRT_PI + lgammafn(b) - lgammafn(.5 + b);
129 xodd = pbeta(x, a, b, /*lower*/TRUE, /*log_p*/FALSE);
130 godd = 2. * rxb * exp(a * log(x) - albeta);
131 tnc = b * x;
132 xeven = (tnc < DBL_EPSILON) ? tnc : 1. - rxb;
133 geven = tnc * rxb;
134 tnc = p * xodd + q * xeven;
135
136 /* repeat until convergence or iteration limit */
137 for(it = 1; it <= itrmax; it++) {
138 a += 1.;
139 xodd -= godd;
140 xeven -= geven;
141 godd *= x * (a + b - 1.) / a;
142 geven *= x * (a + b - .5) / (a + .5);
143 p *= lambda / (2 * it);
144 q *= lambda / (2 * it + 1);
145 tnc += p * xodd + q * xeven;
146 s -= p;
147 /* R 2.4.0 added test for rounding error here. */
148 if(s < -1.e-10) { /* happens e.g. for (t,df,ncp)=(40,10,38.5), after 799 it.*/
149 ML_WARNING(ME_PRECISION, "pnt");
150 #ifdef DEBUG_pnt
151 REprintf("s = %#14.7Lg < 0 !!! ---> non-convergence!!\n", s);
152 #endif
153 goto finis;
154 }
155 if(s <= 0 && it > 1) goto finis;
156 errbd = (double)(2. * s * (xodd - godd));
157 #ifdef DEBUG_pnt
158 REprintf("%3d %#9.4g %#9.4g|%#11.4Lg %#11.4Lg %#11.4Lg %#14.10Lg %#9.4g\n",
159 it, 1e5*(double)godd, 1e5*(double)geven, p, q, s, tnc, errbd);
160 #endif
161 if(fabs(errbd) < errmax) goto finis;/*convergence*/
162 }
163 /* non-convergence:*/
164 ML_WARNING(ME_NOCONV, "pnt");
165 }
166 else { /* x = t = 0 */
167 tnc = 0.;
168 }
169 finis:
170 tnc += pnorm(- del, 0., 1., /*lower*/TRUE, /*log_p*/FALSE);
171
172 lower_tail = lower_tail != negdel; /* xor */
173 if(tnc > 1 - 1e-10 && lower_tail)
174 ML_WARNING(ME_PRECISION, "pnt{final}");
175
176 return R_DT_val(fmin2((double)tnc, 1.) /* Precaution */);
177 }
178