[WP34s] A set of new quantile functions

[Dieter]

Some 34s quantile functions for the statistical distributions do not seem to work as fast and reliably as they could. I recently suggested an update for Student's t quantile, based on a new estimate and the Halley method for fast convergence.

I now would like to propose a similar solution for the Chi² quantile. Here things are a bit more complex: a regular Halley iteration can be set up easily, as shown in another earlier thread, but there are issues with large degrees of freedom. Although the first Chi² estimate is only approx. 11% low in the worst case, this can mean that the CDF of this can be off by several orders of magnitude. This leads to very slow convergence in the first 2...10 iterations which essentially apply the same correction until the estimate is sufficiently accurate and rapid convergence starts.

So I chose a different approach: instead of cdf(x) – p, better solve ln(cdf(x)/p) = 0. This requires a bit of calculus, but it is not too hard to do.

The Student case was already discussed earlier, so the code below is just a slight update. Please note that the incomplete Beta function recently was updated and now requires a new order of its three arguments in X, Y and Z. The code below assumes the previous order.

Compared to the current version included in the 34s, the Normal quantile was streamlined and shortened. It also takes advantage of the new Normal estimate.

The result below is part of a set of programs for the Normal, Student and Chi² quantile. They all rely on a short routine that provides an estimate for the Normal quantile. Compared to the routine previously used in the 34s, only some constants were changed to provide better accuracy. More important, "GNQ" now accepts any probability between 0 and 1 (exclusively) and it returns a correctly signed result. As an additional benefit, min(p, 1–p) is returned in Y which can be handy here and there.

Finally, the following set of programs is what I got until now. This code is to be run in double precision mode, four stack levels are sufficient. As usual, the number of degrees of freedom is assumed in register J. Otherwise the programs only use R00 and R01 and flags 00 and 01. I could not find any severe errors, but please do your own thorough tests and do not consider this perfect. As usual, any error reports and suggestions are welcome.

Code:


*LBL'NQF' // Normal quantile function
XEQ'GNQ'
CF 00
x<0?
SF 00     // save sign in flag 00
x[<->] Y
STO 00    // store min(p, 1-p)
# 002
STO 01    // not more than 2 iterations are required
RCL Z
ABS
FILL
INC X
RSD 03
x[!=]1?
GTO 000
DEC 01    // estimate < 0.005 requires just one iteration
DROP
SDL 016
x>1?
GTO 000
DROP      // estimate < 5E-17 is already exact
GTO 004
*LBL 000
DROP
*LBL 001
FILL
x<1?      // any threshold between 0.254 and 1.28 should be ok
GTO 002
[PHI][sub-u](x)
RCL- 00
GTO 003   // large x: evaluate upper cdf(x) - p
*LBL 002
x[^2]     // small x: evaluate (0.5-p) minus integral from 0 to x
# 1/2     // by means of the incomplete Gamma function
[times]
# 1/2
I[GAMMA][sub-p]
RCL[times] L
+/-
# 1/2
RCL- 00
+
*LBL 003
x[<->] Y
[phi](x)
/
ENTER[^]  // use extrapolation scheme as shown in Abramowitz & Stegun p. 954
x[^3]
RCL T
x[^2]
STO+ X
INC X
[times]
# 006
/
# 1/2
RCL[times] T
RCL[times] Z
RCL[times] Z
+
+
+
DSE 01
GTO 001
*LBL 004
FS?C 00   // adjust sign
+/-       // and exit
END

*LBL'TQF' // Student's t quantile function
ENTER[^]
+/-
INC X
MIN
CF 00
x[!=]? L
SF 00     // save sign in flag 00
STO 00    // store min(p, 1-p)
# 006
STO 01    // do at most 6 iterations (usually 3 or 4)
# 012
RCL J
+/-
y[^x]
RCL 00
x>? Y    // threshold for tail estimate is p = 12^(-dof)
GTO 000
RCL J    // tail estimate for p close to 0 (or 1)
STO+ X
[times]
# [pi]
RCL L
# 1/2
x[^2]
DEC X
+
/
[sqrt]
[times]
RCL J
[^x][sqrt]y
RCL J
[sqrt]
x[<->] Y
/
GTO 001    // start iteration
*LBL 000
XEQ'GNQ'   // estimate for all other p, *not* close to 0 or 1
x[^2]      // t is a tranformation of the Normal z-quantile
# eE
RCL[times] J
1/x
INC X
[times]
RCL/ J
e[^x]-1
RCL[times] J
[sqrt]
*LBL 001
FILL
# 1/2     // use different methods for evaluating cdf(t) - p
x>? Y     // for t close or not close to the center
GTO 002
DROP
t[sub-u](x)
RCL- 00
GTO 003
*LBL 002
DROP
x[^2]     // for t close to the center
ENTER[^]  // evaluate cdf(t) - p by means of the incomplete Beta function
RCL+ J
/
# 1/2
RCL[times] J
RCL L
I[beta]   // IMPORTANT: newer 34s firmware requires different order of arguments! See below.
RCL[times] L
+/-
# 1/2
RCL- 00
+
*LBL 003
RCL T     // use Halley method for fast convergence
t[sub-p](x)
/
ENTER[^]
RCL[times] T
RCL J
INC X
[times]
RCL T
x[^2]
RCL+ J
STO+ X
/
DEC X
/
-
CNVG? 00  // if the last two approximations agree in 14 digits
SKIP 003  // the latest result can be considered converged
DSE 01
GTO 001
ERR 20    // quit with error if no convergence within 6 iterations
FS?C 00
+/-       // else set sign and exit
END

*LBL'CQF' // Chi^2 quantile function
x=0?
RTN       // catch case p = 0
STO 00
# 006
STO 01    // do at most 6 iterations (usually 3...4)
# 019
SDR 001   // threshold for lower tail estimate is 1/pi * 1.9^-n
RCL J
x=1?
DEC X     // handle threshold for 1 dof separately (1/pi is fine)
+/-
y[^x]
# [pi]
/
RCL 00
x<? Y
GTO 000   // jump to lower tail estimate
XEQ'GNQ'  // else get correctly signed Normal quantile
# 222     // and do a Wilson-Hilferty transformation
SDR 003
RCL/ J
STO Z
[sqrt]
[times]
INC X
RCL- Y
x[^3]
RCL[times] J   // Wilson-Hilferty estimate
# eE           // is adjusted for p close to 1
RCL[times] J   // i.e. for Chi^2 > e*dof + 8
# 008
+
x[<->] Y
x<? Y          // Chi^2 not large?
GTO 001        // then start iteration
# 1/2          // otherwise adjust estimate
[times]
LN
# 1/2
RCL[times] J
DEC X
[times]
+/-
RCL 00
+/-
LN1+x
+
# 1/2
RCL[times] J
LN[GAMMA]
+
STO+ X
+/-
GTO 001
*LBL 000   // lower tail estimate
RCL[times] J
# 1/2
[times]
LN
# 1/2
RCL[times] J
LN[GAMMA]
+
STO+ X
RCL/ J
e[^x]
STO+ X
*LBL 001   // iteration starts here
FILL       // uses Halley method for f(x) = ln(cdf(chi^2) / p)
x[>=]? J   // with different methods for low and high Chi^2
GTO 002
[chi][^2]  // low Chi^2
ENTER[^]
RCL/ 00
LN
GTO 003    // provide lower cdf in Y and ln(cdf/p) in X
*LBL 002
# 001
ENTER[^]
RCL- 00
RCL Z
[chi][^2][sub-u]
STO- Z    // Z := 1-Z, i.e. upper cdf => lower cdf
-
RCL/ 00   // provide lower cdf in Y and ln1+x((1-p) - upper_cdf) / p) in X
LN1+x     // the latter is mathematically equivalent to ln(lower_cdf/p)
*LBL 003
x[<->] Y  // start Halley approximation
RCL Z
[chi][^2][sub-p]
x[<->] Y
/
STO Z
/
RCL J
DEC X
DEC X
RCL- T
RCL/ T
[<->] ZXYT  // all this can be done with four stack levels ;-)
STO+ X
-
# 004
/
RCL[times] Y
+/-
INC X
/
-
CNVG? 00  // if the last two approximations agree in 14 digits
SKIP 003  // the latest result can be considered converged
DSE 01
GTO 001
ERR 20    // quit with error if no convergence within 6 iterations
END       // else exit

**LBL'GNQ'  // Provides a guess for the Normal quantile
ENTER[^]    // input: 0 < p < 1
+/-         // output: X = signed quantile, Y = min(p, 1-p)
INC X       // max. absolute and relative error < 0.006
MIN
CF 01
x[!=]? L
SF 01       // save sign in flag 01
# 002
SDR 001     // threshold is 0.2  (was 0.23 before)
x[<->] Y
x[>=]? Y    // use different estimates for p < or > 0.2
GTO 000
FILL        // tail estimate for p < 0.2 or > 0.8
LN
STO+ X
+/-
ENTER[^]
DEC X
# [pi]
[times]
STO+ X
[sqrt]
RCL[times] T
LN
STO+ X
+/-
[sqrt]
x[<->] Y
# 132
STO+ X
SDR 003     // 0.264  (was 0.254 in earlier versions)
x[<->] Y
/
+
GTO 001
*LBL 000
ENTER[^]    // central estimate for 0.2 < p < 0.8
+/-
# 1/2
+
# [pi]
STO+ X
[sqrt]
[times]
ENTER[^]
x[^3]
# 005       // was 6 in earlier versions
/
+
*LBL 001
FS?C 01   // adjust sign
+/-       // negative for p<0.5, positive for p>0.5
END

Edit: there was an update in the incomplete Beta function. If your 34s displays this function as "Iβ" you are running the older version and the original code above works fine. If the function is named I_x instead, this is the newer version with a different order of the three arguments. In this case use the following code between LBL 02 and 03 in "TQF":

Code:

...
*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]
# 1/2
[times]
+/-
# 1/2
RCL- 00
+
*LBL 003
...

And finally: Pauli, what do you think ?-)

Dieter

[walter b]

AFAICS, Pauli uploaded some modifications following your suggestions last night. Will be in next build of WP 34S and WP 31S.

d:-)

[Paul Dale]

These new distribution functions aren't in the image. I'll get to them once things calm down a bit here.

The changes were to the L.R. command.


- Pauli

[Paul Dale]

A new build will depend on Marcus. The L.R. changes are committed.

The distribution changes aren't going to happen immediately. I've more important things to worry about than calculators at present.


- Pauli

[Dieter]

(05-26-2014 09:36 AM)fhub Wrote:  That's a bit strange in my opinion:
First you complain about the illogical order of arguments for IBeta a few weeks ago, then Pauli has changed that (and it was good so), and now you provide a new set of code for distribution functions with exactly the same old (illogical) order of arguments?

Yes. The essential point was not the order of the Iβ arguments, but the fact that the function worked different from what the manual said. This has been corrected. Iβ now works "as advertized", which indeed is more logical.

But there is another point. The beauty of the German language has gifted us with the term "Die normative Kraft des Faktischen". What do you think, how many 34s are currently using the latest firmware? Five percent? Ten? That's why I assumed the previous Iβ implementation. All others may add a simple stack shuffle command.

Quote:It would really be nice to have a (quite) bugfree WP34s version finally

That's why I humbly proposed these three routines. The previous versions sometime threw errors or required inacceptable execution times.

What about you? Your TVM solver has made it into the 34s and added some essential functionality. Great. But what about another project? While I worked on the quantile functions I discovered both the Iβ bug as well as the L.R. problem. Both have been corrected – and the 34s became a bit better. This is a community effort – "ask not what the 34s can do for you – ask what you can do for the 34s project". :-)

Quote:(and I don't care whether any high-level function is accurate to 37 or 'only' 36 digits internally ).

It's more like 30...34. There are limitations in the algorithm and the functions it uses as well as in the hardware. But I think 30 or 32 digits are good enough. If you want all 34 in DP mode, having these functions coded in C is the only way. ;-)

Dieter

[Dieter]

(05-26-2014 10:45 AM)Paul Dale Wrote:  The distribution changes aren't going to happen immediately. I've more important things to worry about than calculators at present.

Take your time. And most important: the suggested code is just a humble proposal. I do not say it is perfect and I strongly suggest you (and the rest of the community) do some thorough tests before this is included. Maybe there is a better approach, a faster solver, a better χ² estimate than good old Wilson-Hilferty, or whatever. I only can say that up to now everything seems to work fine here on my hardware 34s – within the limiations of the respective CDFs and Iβ resp. IΓ. Maybe these could get updated as well...?-)

Dieter

[Dieter]

(05-25-2014 09:40 PM)Dieter Wrote:  Finally, the following set of programs is what I got until now.

Addendum: I updated the original post and added a modified piece of code for use with the newest 34s firmware with the modified incomplete Beta function. I also noticed that this function now appears as I_x, which I think is a good idea.

Dieter

[walter b]

(05-28-2014 01:29 PM)Dieter Wrote:  I updated the original post and added a modified piece of code for use with the newest 34s firmware with the modified incomplete Beta function. I also noticed that this function now appears as I_x, which I think is a good idea.

Thanks.

d:-)

[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

[Dieter]

(05-25-2014 09:40 PM)Dieter Wrote:  I now would like to propose a similar solution for the Chi² quantile. Here things are a bit more complex: a regular Halley iteration can be set up easily, as shown in another earlier thread, but there are issues with large degrees of freedom. Although the first Chi² estimate is only approx. 11% low in the worst case, this can mean that the CDF of this can be off by several orders of magnitude. This leads to very slow convergence in the first 2...10 iterations which essentially apply the same correction until the estimate is sufficiently accurate and rapid convergence starts.

Recent tests showed that a similar behaviour may also occur with the Student quantile function. Up to now I only noticed this in extreme cases out of the SP working range, and here six or seven iterations seem to be sufficient. So the number of allowed iterations initially stored in R01 should be set to 8 for the moment. I will take a closer look at this and see if an approach similar to the Chi² case will eliminate this.

Dieter

[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

[Dieter]

It's me again. Here is the third reply to my own post. #-)

(06-01-2014 01:18 PM)Dieter Wrote:  

(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.
...
If the new code is too long and/or uses too much memory, the iteration limit in R01 may be removed.

Final update: there is another solution. The original (shorter) TQF version can be used with a slight modification.

The program uses two different estimates. One for the central region and another for p close to 0 or 1. At least one of these seems to be good enough to ensure reasonably fast convergence. But in some cases, especially for very high dof values, choosing the one or the other is quite tricky since there is just a narrow corridor of probabilities where both approximations work equally well. So the crucial point is the decision whether to use the central or the tail estimate.

Up to now this was done by checking whether the input is greater or less than 12^–n. This mostly works well, but it is not sufficiently exact for a few critical cases. After some tests I can suggest a different threshold that seems to work fine: simply replace the constant 12 with the term sqrt(n) + 7,5.

So these lines near the beginning of TQF...

Code:

...
# 012
RCL J
+/-
y[^x]
...

...are replaced with these:

Code:

...
RCL J
[sqrt]
# 075
SDR 001
+
RCL J
+/-
y[^x]
...

And that's it.

So there are essentially two ways of solving the problem of slow convergence in some extreme cases:

1. The safe method: use a robust algorithm. This is the approach of the second version I posted, i.e. solving ln(cdf/p) = 0. However, this requires more code and memory.

2. Since the initial estimates are good enough, the threshold for switching between the two may be fine-tuned. With the suggested fix the original shorter version of TQF should work fine. I could not find a case where more than four iterations were required.

Dieter

[Paul Dale]

(06-02-2014 08:32 PM)Dieter Wrote:  It's me again. Here is the third reply to my own post. #-)

I'm reading these, just not acting at present.


- Pauli

[Dieter]

(06-02-2014 10:06 PM)Paul Dale Wrote:  I'm reading these, just not acting at present.

Thank you for your feedback. Today I found an even better way of defining the threshold that works far more accurately than required. Should this be used anyway? Yes, I think so. Why? Because we can. ;-)

So if the updated Student quantile should eventually make it into the 34s, please let me know. But again: take your time.

Dieter