[WP34s] A set of new quantile functions
06-01-2014, 01:18 PM
Post: #11
 Dieter Senior Member Posts: 2,397 Joined: Dec 2013
RE: [WP34s] A set of new quantile functions
(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
 « Next Oldest | Next Newest »

 Messages In This Thread [WP34s] A set of new quantile functions - Dieter - 05-25-2014, 09:40 PM RE: [WP34s] A set of new quantile functions - walter b - 05-26-2014, 07:40 AM RE: [WP34s] A set of new quantile functions - Paul Dale - 05-26-2014, 08:06 AM RE: [WP34s] A set of new quantile functions - Paul Dale - 05-26-2014, 10:45 AM RE: [WP34s] A set of new quantile functions - Dieter - 05-26-2014, 06:52 PM RE: [WP34s] A set of new quantile functions - Dieter - 05-29-2014, 12:43 PM RE: [WP34s] A set of new quantile functions - Dieter - 05-28-2014, 01:29 PM RE: [WP34s] A set of new quantile functions - walter b - 05-28-2014, 02:19 PM RE: [WP34s] A set of new quantile functions - Dieter - 05-30-2014, 06:13 AM RE: [WP34s] A set of new quantile functions - Dieter - 06-01-2014 01:18 PM RE: [WP34s] A set of new quantile functions - Dieter - 06-02-2014, 08:32 PM RE: [WP34s] A set of new quantile functions - Paul Dale - 06-02-2014, 10:06 PM RE: [WP34s] A set of new quantile functions - Dieter - 06-03-2014, 08:35 PM

User(s) browsing this thread: 1 Guest(s)