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