[WP34s] Regularized incomplete Beta function
05-04-2014, 05:26 PM (This post was last modified: 05-04-2014 05:33 PM by Dieter.)
Post: #21
 Dieter Senior Member Posts: 2,397 Joined: Dec 2013
RE: [WP34s] Regularized incomplete Beta function
(05-04-2014 06:29 AM)Paul Dale Wrote:  I can only agree with you here. I'm using a fairly generic Newton based solver for all the quantile functions except the normal QF. I know it can be very slow, that is the price for having a single general solver.

Even a simple Newton solver should converge quadratically. I think the slow execution of the current 34s quantile implementations is due to an error in the initial guess for the solver. It should be within approx. 10% of the true quantile. So even a Newton solver should not require much more than five iterations, which means 6 - 10 s in user code.

Quote:The 34S certainly doesn't have space for a C implementation of these. I originally implemented all the distribution functions in C but had to rewrite them as keystroke programs for space reasons. Adding additional custom solvers is likewise going to consume precious bytes. Thus, the 34S is unlikely to see any distribution speed ups. However, the 31S quite possibly will.

Great. And many things are easier to do as there is no DP mode. ;-)

Quote:I'd like to rewrite the distribution code there in native C -- I can't just drop in the original distribution code, we moved on algorithmically.
So, yes I'd love to hear your thoughts on improving these functions.

I recently tried some improvements, especially with the solver. The following ideas refer to the Student distribution, but they should be applicable for the Chi² and Fisher cases as well. The following thoughts assume the quantile is ≥ 0 and p is the upper tail probability (i.e. p ≤ 0.5). As usual, the rest is done with a simple transformation, e.g. QF(1–p) = –QF(p).

So, what can be done?

1. Use a Halley solver. It converges very fast and it can be implemented easily since the second derivative is a function of the first (i.e. the pdf).

Student:
f"(x) = -(n+1)x / (n+x^2) * pdf

Chi²:
f"(x) = (n-x-2)/(2x) * pdf

Fisher:
f"(x) = -(n(m(x-1)+2) + 2mx) / (2x(mx+n)) * pdf

where n (resp. m and n) are the degrees of freedom.

So f"(x) simply is the pdf times a factor. This way a Halley iteration e.g. for the Student case can be written as follows:

r := (StudentUpperCDF(x) - p) / StudentPDF
d := r / ( r * (n+1)x / 2(n+x^2) - 1 )
x := x - d

Due to the roughly cubic convergence, the iteration usually may quit as soon as |d| < 1E–10*x for 30+ digit accuracy. In my 34s program I tried a conservative CNVG 00 (i.e. rel. error < 1E–14) which usually returns results that are as good as a user code program gets when running in DP mode (approx 30-34 digits).

The same idea can be used with the Chi² and Fisher quantile.

2. I slightly modified the initial guess for the Chi² quantile and I tried a new approach in the Student case, based on a 1992 paper on bounds of various quantiles by Fujikoshi and Mukaihata. The idea is a simple transformation of the normal quantile z:

t = sqrt(n * (ea * z^2/n - 1 ))

Where a is close to 1 and a function of n (or simply 1 so that it can be omitted). I used a = 1 + 1/(e*n) which works very well with the Normal estimate I used (simply the one in the Normal quantile function).

This works fine for the center but less so for the distribution tails. For all p < 12–n I use a slight modification of the tail approximation suggested years ago:

u = 2 * p * n * sqrt(pi / (2 * n - 0.75))
t = sqrt(n) / u ^ (1 / n)

The 0.75 originally was a 1, but changing this value improves the results for low n.

Although the Student estimate originally was intended for n≥3 it also works well for n=2 or even n=1. Usually it converges within three iterations, here and there in maybe four. This means that the code for n=1 and n=2 (direct solutions) may be omitted.

For x close to 0 (i.e. p near 0.5) the expression t_u(x) – p loses accuracy due to digit cancellation. So I used the same idea as in the Normal quantile routine and had this value calculated differently for small t, using the incomplete Beta function. Yes, that's why I found the bug discussed in this thread. ;-)

Code:
 LBL 'TQF'  ' T Quantile Function ENTER      ' probability in X, degrees of freedom in register J +/- INC X MIN CF 00 x!=? L SF 00       ' set flag 00 if p > 0.5 CF 01       ' clear error flag STO 00      ' save p in R00 # 005 STO 01      ' do not more than 5 iterations # 012 RCL J +/- y^x RCL 00 x>? Y GTO 00 RCL J       ' estimate for small p STO+ X × # pi RCL L # 1/2 x² RCL+ L    ' = 0.75 - / sqrt × RCL J xrooty RCL J sqrt x<> Y / GTO 01 LBL 00     ' estimate for low and moderate t XEQ 'GNQ'  ' get guess for the normal quantile x² # eE RCL× J 1/x INC X × RCL/ J e^x-1 RCL× J sqrt LBL 01     ' iteration starts here FILL # 1/2 x>? Y GTO 02 DROP t_u(x) RCL- 00 GTO 03 LBL 02 DROP x² ENTER RCL+ J / # 1/2 RCL× J RCL L IBeta RCL× L +/- # 1/2 RCL- 00 +          ' cdf(t) - p = (0,5-p) - 1/2 IBeta(x=t^2/(n+t^2), a=1/2, b=n/2) LBL 03 RCL T t_p(x) / ENTER RCL× T RCL J INC X × RCL T x² RCL+ J STO+ X / DEC X / - CNVG? 00 SKIP 003 DSE 01 GTO 01 SF 01       ' Raise error flag if no convergence after 5 iterations FS?C 00 +/-         ' adjust sign FS?C 01 ERR 20      ' if error, display "no root found" and exit with last approximation END LBL 'GNQ'   ' input: p =< 0.5 # 232       ' output: Normal estimate > 0 SDR 003 x<>Y x>? Y GTO 00 FILL        'Normal estimate for p up to 0.232 LN STO+ X +/- ENTER DEC X # pi × STO+ X sqrt RCL× T LN STO+ X +/- sqrt x<>Y # 004 × 1/x + RTN LBL 00   ' Normal estimate for p close to the center +/- # 1/2 + # pi STO+ X sqrt × ENTER x³ # 006 / + RTN END

For best accuracy this should run in DP mode. The program exits if the last two approximations agree in approx. 14 digits. At this point the result usually carries 30+ valid digits.

Here are some examples:

Code:
 10 STO J    0,1 XEQ"TQF" => -1,372183641110335627219156967662554 in 3,2 s    exact result:   -1,37218364111033562721915696766255392  1E-20 XEQ"TQF" => -256,4346993185261855315362349874343 in 3,1 s    exact result:   -256.434699318526185531536234987434334  0,5 ENTER 1E-16 -        XEQ"TQF" => -2,569978034930492409497513483729480 E-16 in 1,9 s    exact result:   -2,56997803493049240949751348372947856 E-16  1 STO J  0,05  XEQ"TQF" => -6,313751514675043098979464244768186 in 3 s    exact result:   -6,3137515146750430989794642447681860594  1E-10 XEQ"TQF" => -3183098861,837906715272955512330630 in 1,9 s    exact result:   -3183098861.837906715272955512330627466 100 STO J  0,025 XEQ"TQF" => -1,983971518523552286595184867990389 in 5,1 s    exact result:   -1,983971518523552286595184867990339165

Please note that the results for n=1 are exact to 33 resp. 34 digits while the current implementation (that calculates the result directly) gets only 32 resp. 24 (!) digits right. Either the internal tangent function is not that accurate or the current implementation does not evaluate the quantile as 1/tan(1E–10*180°) which would yield a nearly perfect result. ;-)

Dieter
 « Next Oldest | Next Newest »

 Messages In This Thread [WP34s] Regularized incomplete Beta function - Dieter - 05-01-2014, 07:17 PM RE: [WP34s] Regularized incomplete Beta function - Thomas Klemm - 05-01-2014, 07:48 PM RE: [WP34s] Regularized incomplete Beta function - Dieter - 05-01-2014, 08:00 PM RE: [WP34s] Regularized incomplete Beta function - Thomas Klemm - 05-01-2014, 08:28 PM RE: [WP34s] Regularized incomplete Beta function - walter b - 05-01-2014, 09:08 PM RE: [WP34s] Regularized incomplete Beta function - walter b - 05-01-2014, 09:23 PM RE: [WP34s] Regularized incomplete Beta function - Dieter - 05-01-2014, 09:26 PM RE: [WP34s] Regularized incomplete Beta function - Paul Dale - 05-01-2014, 09:55 PM RE: [WP34s] Regularized incomplete Beta function - walter b - 05-01-2014, 10:14 PM RE: [WP34s] Regularized incomplete Beta function - Thomas Klemm - 05-01-2014, 11:13 PM RE: [WP34s] Regularized incomplete Beta function - Thomas Klemm - 05-01-2014, 11:01 PM RE: [WP34s] Regularized incomplete Beta function - Dieter - 05-02-2014, 08:48 PM RE: [WP34s] Regularized incomplete Beta function - walter b - 05-02-2014, 09:25 PM RE: [WP34s] Regularized incomplete Beta function - Thomas Klemm - 05-02-2014, 10:22 PM RE: [WP34s] Regularized incomplete Beta function - Dieter - 05-02-2014, 10:34 PM RE: [WP34s] Regularized incomplete Beta function - walter b - 05-02-2014, 10:46 PM RE: [WP34s] Regularized incomplete Beta function - Thomas Klemm - 05-02-2014, 11:29 PM RE: [WP34s] Regularized incomplete Beta function - Dieter - 05-02-2014, 09:03 PM RE: [WP34s] Regularized incomplete Beta function - Paul Dale - 05-04-2014, 06:29 AM RE: [WP34s] Regularized incomplete Beta function - Dieter - 05-04-2014 05:26 PM RE: [WP34s] Regularized incomplete Beta function - Paul Dale - 05-05-2014, 02:20 AM RE: [WP34s] Regularized incomplete Beta function - Dieter - 05-05-2014, 09:13 PM RE: [WP34s] Regularized incomplete Beta function - Paul Dale - 05-05-2014, 09:57 PM RE: [WP34s] Regularized incomplete Beta function - Paul Dale - 05-05-2014, 08:26 AM RE: [WP34s] Regularized incomplete Beta function - Dieter - 05-11-2014, 01:50 PM RE: [WP34s] Regularized incomplete Beta function - Paul Dale - 05-11-2014, 09:51 PM RE: [WP34s] Regularized incomplete Beta function - Dieter - 05-12-2014, 06:18 AM RE: [WP34s] Regularized incomplete Beta function - Paul Dale - 05-12-2014, 06:58 AM RE: [WP34s] Regularized incomplete Beta function - Dieter - 05-12-2014, 07:47 PM RE: [WP34s] Regularized incomplete Beta function - walter b - 05-12-2014, 08:14 PM RE: [WP34s] Regularized incomplete Beta function - Paul Dale - 05-12-2014, 09:59 PM RE: [WP34s] Regularized incomplete Beta function - Dieter - 05-22-2014, 08:51 PM RE: [WP34s] Regularized incomplete Beta function - Dieter - 05-25-2014, 09:43 PM RE: [WP34s] Regularized incomplete Beta function - Manolo Sobrino - 05-03-2014, 01:38 PM RE: [WP34s] Regularized incomplete Beta function - Paul Dale - 05-05-2014, 11:02 AM RE: [WP34s] Regularized incomplete Beta function - Dieter - 05-05-2014, 08:30 PM RE: [WP34s] Regularized incomplete Beta function - Paul Dale - 05-05-2014, 09:49 PM RE: [WP34s] Regularized incomplete Beta function - Dieter - 05-05-2014, 10:04 PM RE: [WP34s] Regularized incomplete Beta function - Paul Dale - 05-05-2014, 10:16 PM RE: [WP34s] Regularized incomplete Beta function - Dieter - 05-05-2014, 11:09 PM RE: [WP34s] Regularized incomplete Beta function - Paul Dale - 05-06-2014, 12:54 AM RE: [WP34s] Regularized incomplete Beta function - Dieter - 05-06-2014, 11:01 PM RE: [WP34s] Regularized incomplete Beta function - Paul Dale - 05-07-2014, 12:12 AM RE: [WP34s] Regularized incomplete Beta function - Dieter - 05-08-2014, 06:26 PM RE: [WP34s] Regularized incomplete Beta function - Paul Dale - 05-08-2014, 10:28 PM RE: [WP34s] Regularized incomplete Beta function - Dieter - 05-08-2014, 10:58 PM RE: [WP34s] Regularized incomplete Beta function - Paul Dale - 05-08-2014, 11:29 PM RE: [WP34s] Regularized incomplete Beta function - Dieter - 05-09-2014, 09:06 PM RE: [WP34s] Regularized incomplete Beta function - Paul Dale - 05-10-2014, 12:03 AM RE: [WP34s] Regularized incomplete Beta function - Dieter - 05-11-2014, 01:23 PM RE: [WP34s] Regularized incomplete Beta function - Paul Dale - 05-11-2014, 10:33 PM RE: [WP34s] Regularized incomplete Beta function - Paul Dale - 05-05-2014, 11:15 AM RE: [WP34s] Regularized incomplete Beta function - Dieter - 05-05-2014, 08:14 PM RE: [WP34s] Regularized incomplete Beta function - Paul Dale - 05-05-2014, 09:40 PM RE: [WP34s] Regularized incomplete Beta function - Paul Dale - 05-12-2014, 06:05 AM

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