Re: A first estimate for Student's t quantile (not only) on the WP34s Message #23 Posted by Oliver Unter Ecker on 3 May 2011, 6:54 a.m., in response to message #22 by Dieter
Thank you very much, Dieter!
Here's my JavaScript "production" code, based on yours:
"StudentTQuantile": function(df, x) {
if (df != Math.floor(df)) throw Error("bad arg");
if (x <= 0) return +Infinity; else if (x >= 1) return -Infinity; else if (x == 0.5) return 0;
var t, p = Math.min(x, 1-x), s = this.sign(0.5 - x);
if (-this.ln(p) < 1.7 * df) {
if (p < 0.2) {
var u = -2*this.ln(p);
x = Math.sqrt(-2*this.ln(p*Math.sqrt((u-1)*2*Math.PI))) + 1/(u*u);
}
else {
var u = (0.5-p)*Math.sqrt(2*Math.PI);
x = u + (u*u*u)/6;
}
var x3 = x*x*x;
u = (x3 + x) / 4 / df;
v = (x3*x*x / 12 + x3 / 4) / (df*df);
t = x + u + v;
}
else {
u = 2*p*df*Math.sqrt(Math.PI/(2*df - 1));
t = Math.sqrt(df) / Math.pow(u,(1/df));
}
var d, k = 10;
do {
d = (this.StudentTDistribution(df, t) - p) / this.StudentTPDF(df, t);
t += d;
}
while (d*d > 1e-30 && --k);
return t*s;
}
You'll notice the change in the termination criterion to get the result within one ULP (for IEEE754 doubles) but providing a hard iteration max to deal with too slow convergence and oscillations.
Here the results for your extreme values test:
p n exact result ND1 (max 10 iters) ND1 (max 100 iters)
------------------------------------------------------------------------------------
1 E-5 10 7.526954019339 7.526954019333 7.526954019334
1 E-10 10 25.466008 25.466010 25.466008
1 E-15 10 81.040891 80.3924 81.05831
1 E-20 10 256.434699 -Infinity -Infinity
1 E-40 10 2.564 E+05 -Infinity -Infinity
1 E-100 10 2.564 E+10 19881486106.20 -Infinity
I guess a PDF at the mercy of the gamma function, may, in part, explain the slow convergence.
I'm going with the max of 10 iterations. For non-extreme values, results within +/- 1 ULP seem to be found in typically 3-5 iterations. The first estimate is great!
Do you have an iOS device? I'd love to give you a redeem code for a copy of ND1 (and CalcPad when it comes out).
Cheers.
|