[WP34s] A set of new quantile functions

[Dieter]

(05-26-2014 10:45 AM)Paul Dale Wrote:  The distribution changes aren't going to happen immediately.

That's fine since another update may be required as well.

I did some further tests with the recently suggested Student quantile function and found that close to zero the CDF of the calculated quantile often significantly differed from the original probability. Some tests with 42-digit precision showed that the quantile was fine (33 or 34 correct digits), while there are some issues with the WP34s t(x) function.

Near t=0 the Student CDF should be roughly 0,5 + r·t where r is a factor between 1/pi and 1/sqrt(2·pi), i.e. somewhere between 0,3 and 0,4. So the double precision results for the Student CDF should look like this:

Code:

    t          CDF(t)
-----------------------------------------------------------------
   0,1         0,53...
   0,01        0,503...
   0,001       0,5003...
   1E-5        0,50000 3...
   1E-10       0,50000 00000 3...
   1E-20       0,50000 00000 00000 00000 3...
   1E-33       0,50000 00000 00000 00000 00000 00000 0003 or ...4
=< 1E-34       0,50000 00000 00000 00000 00000 00000 0000
-----------------------------------------------------------------

But actually there is a much larger t below which a plain 0,5 is returned:

Code:

    t          WP34s t(x) function for 10 dof
-----------------------------------------------------------------
   1E-16       0,50000 00000 00000 03891 08383 96603 1050
   1E-17       0,50000 00000 00000 00000 00000 00000 0000
-----------------------------------------------------------------

For further investigations I used a little program that provides a 34-digit t-value of the same magnitude as entered (pi/3 * input):

Code:

LBL A
 pi
 x
 3
 /
ENTER
 x²
ENTER
RCL+J
 /
 1/2
ENTER
RCLxJ
x<> Z
 Ix     // with recent firmware
 1/2
 x
RCL+L  // = exact CDF for t close to zero
x<> Y
t(x)   // = CDF as given by the current 34s function
 -     // return absolute error
RTN

10 STO J  // set 10 degrees of freedom

0,1    [A]   1,1 E-33
0,01   [A]  -5,81 E-32
0,001  [A]  -4,085 E-31
1E-5   [A]  -4,2508 E-29
1E-10  [A]   2,8046 E-24
1E-15  [A]  -6,2697 E-19
1E-16  [A]   1,8365 E-18
5E-17  [A]   2,0374 E-17
1E-17  [A]   4,0747 E-18
1E-18  [A]   4,0747 E-19
1E-19  [A]   4,0747 E-20
...
1E-29  [A]   4,0747 E-30
1E-30  [A]   4,075 E-31
1E-31  [A]   4,07 E-32
1E-32  [A]   4,1 E-33
1E-33  [A]   4 E-34
1E-34  [A]   0

Since the true result is approx. 0,5 the returned error can be read as "units in the xth digit". So for t near 5 E–17 the error is 2 units in the 17th digit. In other words, accuracy degrades down to single precision level. I assume this happens due to digit cancellation when the argument x = n/(t²+n) for the incomplete Beta function is evaluated, since for very small t this rounds to 1. That's why I suggest that for small t the Student CDF should be evaluated using 1–x = t²/(t²+n), as shown in the routine above.

BTW the recently posted Student quantile function is not affected by all this as it already uses the suggested method for |t| < 0,5.

EDIT:

I looked at the current Student CDF code in t.wp34s, and indeed it works as assumed. Here is a fix that uses a different method for small t less than 1:

Code:

                        ...
cdf_t_return::          x<1?
                        JMP cdf_t_small
                        RCL J
                        SWAP 
                        x[^2]
                        RCL+ J
                        /
                        Num 1/2
                        RCL[times] J
                        Num 1/2
                        x[<->] Z
                        I[sub-x]
                        Num 1/2
                        [times]
                        RTN

cdf_t_small::           x[^2]
                        ENTER[^]
                        RCL+ J
                        /
                        Num 1/2
                        ENTER[^]
                        RCL[times] J
                        x[<->] Z
                        I[sub-x]
                        Num 1/2
                        [times]
                        +/-
                        RCL+ L
                        RTN

cdf_t_invert::          GSB cdf_t_return
                        +/-  
                        INC X
                        RTN

Dieter