[WP34s] Regularized incomplete Beta function

05042014, 05:26 PM
(This post was last modified: 05042014 05:33 PM by Dieter.)
Post: #21




RE: [WP34s] Regularized incomplete Beta function
(05042014 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. 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) = (nx2)/(2x) * pdf Fisher: f"(x) = (n(m(x1)+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 3034 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 * (e^{a * 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:
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:
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 »

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