[WP34s] A set of new quantile functions

[Dieter]

(05-30-2014 06:13 AM)Dieter Wrote:  Recent tests showed that a similar behaviour may also occur with the Student quantile function.
(...)
I will take a closer look at this and see if an approach similar to the Chi² case will eliminate this.

It does. Here is an improved version of the Student quantile function that uses the same approach as the Chi² quantile: it solves ln(cdf/p) = 0. It does not require more than four stack levels. With eight levels the whole thing might get shortened a bit.

The 34s code is rather complex and hard to understand, so here it is in a more readable form:

Code:


imax = 5
q = min(p, 1-p)
t = guess_t(q, n)

For i = 1 To imax
    
    u = t / (n + t * t)
    r = (n + 1) * u
    pdf = tpdf(t, n)
       
    If t >= 0.325 Then  ' here the cdf is < 0.4 for n = 1...\infty
       cdf = tcdf(t, n) ' upper tail Student CDF
       s = pdf / cdf
       r = r - s
       a = Ln(cdf / q) / s
    Else
       cdf = 0.5 * ((1 - 2*q) - betai(t*u, 1/2, n/2))
       a = cdf / pdf
    End If
    
    dt = a / (a * r / 2 - 1)
    t = t - dt

Next

return t * sign(p-0.5)

The 34s code is the same as before, only the iteration loop at LBL 01 is replaced with this:

Code:

...
*LBL 001
FILL
# 1/2
x>? Y           // 0.325 would be perfect, but 0.5 or even 1 will do either
GTO 002
SIGN
RCL+ J          // n+1
[times]
x[<->] Y
x[^2]
RCL+ J
/               // r := f"/f' = (n+1)t²/(n+t²)
x[<->] Y
t[sub-u](x)     // cdf = upper tail CDF
RCL L
t[sub-p](x)
RCL/ Y          // pdf/cdf   cdf  r  t
x[<->] Y
RCL/ 00
LN              // ln(cdf/p)  pdf/cdf  r  t
x[<->] Y
/               // -f/f' = ln(cdf/p) * cdf/pdf
x[<->] Y        // r  ln(cdf/p)*cdf/pdf  t  t
RCL- L          // r - pdf/cdf
# 1/2
[times]         // (r - pdf/cdf)/2 = f"/f'/2
RCL[times] Y
GTO 003

*LBL 002
RCL[times] J    // X was 1/2, so this is n/2
x[<->] Y
x[^2]
ENTER[^]
RCL+ J
/               // x = t²/(t²+n)
# 1/2
x[<->] Y        // stack now is  x  1/2  n/2  t
I[sub-x]
+/-
RCL 00
STO+ X
+/-
INC X
+
# 1/2
[times]       // = f(x) = upper tail cdf(x) - p
x[<->] Y
t[sub-p](x)
/             // = -f(x)/f'(x)
ENTER
RCL[times] T  // evaluate -f/f' * f"/f'/2
RCL J         // = -f/f' * (n+1)t/2(t²+n)
INC X
[times]
RCL T
x[^2]
RCL+ J
STO+ X
/

*LBL 003
DEC X
/
-
CNVG? 00
SKIP 003
DSE 01
GTO 001
ERR 20
FS?C 00
+/-
RTN

As before, 5 iterations are initially stored in R01. If the new code is too long and/or uses too much memory, the iteration limit in R01 may be removed. More than 5 iterations are only required in a limited domain for very small p and high dof. Up to now I have not found any cases within the SP working range. So if some more iterations are OK in these cases, the original code may be used as well. The worst I've seen so far are 12 loops (at 5000 dof, p=10^-5400). However, the new code does it in three. ;-)

Dieter