1 module nt;
2
3 import rmathd.common;
4 import rmathd.t;
5 import rmathd.normal;
6 import rmathd.beta;
7 import rmathd.chisq;
8
9
10 /*
11 ** toms708.d
12 ** dmd nt.d common.d t.d normal.d beta.d chisq.d && ./nt
13 */
14
15
16 T dnt(T)(T x, T df, T ncp, int give_log)
17 {
18 T u;
19 mixin R_D__0!give_log;
20
21 if (isNaN(x) || isNaN(df))
22 return x + df;
23
24
25 /* If non-positive df then error */
26 if (df <= 0.0)
27 return T.nan;
28
29 if(ncp == 0.0)
30 return dt!T(x, df, give_log);
31
32 /* If x is infinite then return 0 */
33 if(!isFinite(x))
34 return R_D__0;
35
36 /* If infinite df then the density is identical to a
37 normal distribution with mean = ncp. However, the formula
38 loses a lot of accuracy around df=1e9
39 */
40 if(!isFinite(df) || df > 1e8)
41 return dnorm!T(x, ncp, 1., give_log);
42
43 /* Do calculations on log scale to stabilize */
44
45 /* Consider two cases: x ~= 0 or not */
46 if (fabs(x) > sqrt(df * DBL_EPSILON)) {
47 u = log(df) - log(fabs(x)) +
48 log(fabs(pnt!T(x*sqrt((df + 2)/df), df + 2, ncp, 1, 0) -
49 pnt!T(x, df, ncp, 1, 0)));
50 /* FIXME: the above still suffers from cancellation (but not horribly) */
51 }
52 else { /* x ~= 0 : -> same value as for x = 0 */
53 u = lgammafn!T((df+1)/2) - lgammafn!T(df/2)
54 - (M_LN_SQRT_PI + .5*(log(df) + ncp*ncp));
55 }
56
57 return (give_log ? u : exp(u));
58 }
59
60
61
62 T pnt(T)(T t, T df, T ncp, int lower_tail, int log_p)
63 {
64 T albeta, a, b, del, errbd, lambda, rxb, tt, x;
65 real geven, godd, p, q, s, tnc, xeven, xodd;
66 int it, negdel;
67
68 /* note - itrmax and errmax may be changed to suit one's needs. */
69
70 const int itrmax = 1000;
71 const static T errmax = 1.0e-12;
72
73 if (df <= 0.0)
74 return T.nan;
75 if(ncp == 0.0)
76 return pt!T(t, df, lower_tail, log_p);
77
78 if(!isFinite(t))
79 return (t < 0) ? R_DT_0!T(lower_tail, log_p) : R_DT_1!T(lower_tail, log_p);
80 if (t >= 0.) {
81 negdel = 0; tt = t; del = ncp;
82 } else {
83 /* We deal quickly with left tail if extreme,
84 since pt(q, df, ncp) <= pt(0, df, ncp) = \Phi(-ncp) */
85 if (ncp > 40 && (!log_p || !lower_tail))
86 return R_DT_0!T(lower_tail, log_p);
87 negdel = 1; tt = -t; del = -ncp;
88 }
89
90 if (df > 4e5 || del*del > 2*M_LN2*(-(DBL_MIN_EXP))) {
91 /*-- 2nd part: if del > 37.62, then p=0 below
92 FIXME: test should depend on `df', `tt' AND `del' ! */
93 /* Approx. from Abramowitz & Stegun 26.7.10 (p.949) */
94 s = 1./(4.*df);
95
96 return pnorm!T(cast(T)(tt*(1. - s)), del,
97 sqrt(cast(T) (1. + tt*tt*2.*s)),
98 lower_tail != negdel, log_p);
99 }
100
101 /* initialize twin series */
102 /* Guenther, J. (1978). Statist. Computn. Simuln. vol.6, 199. */
103
104 x = t * t;
105 rxb = df/(x + df);/* := (1 - x) {x below} -- but more accurately */
106 x = x / (x + df);/* in [0,1) */
107 //#ifdef DEBUG_pnt
108 // REprintf("pnt(t=%7g, df=%7g, ncp=%7g) ==> x= %10g:",t,df,ncp, x);
109 //#endif
110 if (x > 0.) {/* <==> t != 0 */
111 lambda = del * del;
112 p = .5 * exp(-.5 * lambda);
113 //#ifdef DEBUG_pnt
114 // REprintf("\t p=%10Lg\n",p);
115 //#endif
116 if(p == 0.) { /* underflow! */
117
118 /*========== really use an other algorithm for this case !!! */
119 //ML_ERROR(ME_UNDERFLOW, "pnt");
120 //ML_ERROR(ME_RANGE, "pnt"); /* |ncp| too large */
121 //assert(0, "underflow error in pnt");
122 //assert(0, "range error in pnt");
123 return R_DT_0!T(lower_tail, log_p);
124 }
125 //#ifdef DEBUG_pnt
126 // REprintf("it 1e5*(godd, geven)| p q s"
127 // /* 1.3 1..4..7.9 1..4..7.9|1..4..7.901 1..4..7.901 1..4..7.901 */
128 // " pnt(*) errbd\n");
129 // /* 1..4..7..0..34 1..4..7.9*/
130 //#endif
131 q = M_SQRT_2dPI * p * del;
132 s = .5 - p;
133 /* s = 0.5 - p = 0.5*(1 - exp(-.5 L)) = -0.5*expm1(-.5 L)) */
134 if(s < 1e-7)
135 s = -0.5 * expm1(-0.5 * lambda);
136 a = .5;
137 b = .5 * df;
138 /* rxb = (1 - x) ^ b [ ~= 1 - b*x for tiny x --> see 'xeven' below]
139 * where '(1 - x)' =: rxb {accurately!} above */
140 rxb = pow(rxb, b);
141 albeta = M_LN_SQRT_PI + lgammafn!T(b) - lgammafn!T(.5 + b);
142 xodd = pbeta!T(x, a, b, /*lower*/1, /*log_p*/0);
143 godd = 2. * rxb * exp(a * log(x) - albeta);
144 tnc = b * x;
145 xeven = (tnc < DBL_EPSILON) ? tnc : 1. - rxb;
146 geven = tnc * rxb;
147 tnc = p * xodd + q * xeven;
148
149 /* repeat until convergence or iteration limit */
150 for(it = 1; it <= itrmax; it++) {
151 a += 1.;
152 xodd -= godd;
153 xeven -= geven;
154 godd *= x * (a + b - 1.) / a;
155 geven *= x * (a + b - .5) / (a + .5);
156 p *= lambda / (2 * it);
157 q *= lambda / (2 * it + 1);
158 tnc += p * xodd + q * xeven;
159 s -= p;
160 /* R 2.4.0 added test for rounding error here. */
161 if(s < -1.0e-10) { /* happens e.g. for (t,df,ncp)=(40,10,38.5), after 799 it.*/
162 //ML_ERROR(ME_PRECISION, "pnt");
163 //assert(0, "precision error int pnt");
164
165 //#ifdef DEBUG_pnt
166 // REprintf("s = %#14.7Lg < 0 !!! ---> non-convergence!!\n", s);
167 //#endif
168 goto finis;
169 }
170 if(s <= 0 && it > 1) goto finis;
171 errbd = cast(T)(2. * s * (xodd - godd));
172 //#ifdef DEBUG_pnt
173 // REprintf("%3d %#9.4g %#9.4g|%#11.4Lg %#11.4Lg %#11.4Lg %#14.10Lg %#9.4g\n",
174 // it, 1e5*(double)godd, 1e5*(double)geven, p, q, s, tnc, errbd);
175 //#endif
176 if(fabs(errbd) < errmax) goto finis;/*convergence*/
177 }
178 /* non-convergence:*/
179 //ML_ERROR(ME_NOCONV, "pnt");
180 assert(0, "non-convergence in pnt");
181 } else { /* x = t = 0 */
182 tnc = 0.;
183 }
184 finis:
185 tnc += pnorm!T(- del, 0., 1., /*lower*/1, /*log_p*/0);
186
187 lower_tail = lower_tail != negdel; /* xor */
188 if(tnc > 1 - 1e-10 && lower_tail){
189 //ML_ERROR(ME_PRECISION, "pnt{final}");
190 assert(0, "precision error in pnt");
191 }
192
193 return R_DT_val!T(fmin2!T(cast(T)tnc, 1), lower_tail, log_p);
194 }
195
196 /* Currently returns the wrong answer! */
197 T qnt(T)(T p, T df, T ncp, int lower_tail, int log_p)
198 {
199 const static T accu = 1e-13;
200 const static T Eps = 1e-11; /* must be > accu */
201
202 T ux, lx, nx, pp;
203
204 if (isNaN(p) || isNaN(df) || isNaN(ncp))
205 return p + df + ncp;
206
207 /* Was
208 * df = floor(df + 0.5);
209 * if (df < 1 || ncp < 0) ML_ERR_return_NAN;
210 */
211 if (df <= 0.0)
212 return T.nan;
213
214 if(ncp == 0.0 && df >= 1.0)
215 return qt!T(p, df, lower_tail, log_p);
216
217 immutable(T) NINF = -T.infinity;
218 immutable(T) INF = T.infinity;
219 mixin (R_Q_P01_boundaries!(p, NINF, INF)());
220
221 if (!isFinite(df)) // df = Inf ==> limit N(ncp,1)
222 return qnorm!T(p, ncp, 1., lower_tail, log_p);
223
224 mixin R_DT_qIv!p;
225 p = R_DT_qIv;
226
227 /* Invert pnt(.) :
228 * 1. finding an upper and lower bound */
229 if(p > 1 - DBL_EPSILON)
230 return INF;
231 pp = fmin2!T(1 - DBL_EPSILON, p * (1 + Eps));
232 for(ux = fmax2!T(1., ncp);
233 ux < DBL_MAX && pnt!T(ux, df, ncp, 1, 0) < pp;
234 ux *= 2){}
235 pp = p * (1 - Eps);
236 for(lx = fmin2(-1., -ncp);
237 (lx > -DBL_MAX) && (pnt!T(lx, df, ncp, 1, 0) > pp);
238 lx *= 2){}
239
240 /* 2. interval (lx,ux) halving : */
241 do {
242 nx = 0.5 * (lx + ux); // could be zero
243 if (pnt!T(nx, df, ncp, 1, 0) > p)
244 ux = nx;
245 else
246 lx = nx;
247 }
248 while ((ux - lx) > accu * fmax2!T(fabs(lx), fabs(ux)));
249
250 return 0.5 * (lx + ux);
251 }
252
253
254 T rnt(T)(T df, T ncp)
255 {
256 return rnorm!T(ncp, 1.)/sqrt(rchisq!T(df)/df);
257 }
258
259
260 void test_nt()
261 {
262 import std.stdio: writeln;
263 writeln("dnt: ", dnt(2., 40., 2., 0));
264 writeln("pnt: ", pnt(3., 40., 2., 1, 0));
265 writeln("pnt (wrong answer): ", pnt(-3., 40., 2., 1, 0));
266 //writeln("qnt: ", qnt(0.7, 40., 1., 1, 0));
267 writeln("rnt: ", rnt(40., 2.), ", rnt: ", rnt(40., 2.), ", rnt: ", rnt(40., 2.));
268 }
269
270
271