[WP34s] Regularized incomplete Beta function

[Dieter]

(05-06-2014 12:54 AM)Paul Dale Wrote:  Wouldn't we have to force one additional iteration after achieving the convergence precision to guarantee this level of accuracy? A value that differed from the previous by, say, 1E-15, would pass the convergence test yet not be good.

If the last correction was 1E-15, this was the error of the previous approximation. Since Halley's method converges roughly cubic, the expected remaining error after this correction is below 1E-40.

Quote:I'm thinking along the lines of:

Code:

FS? .02
JMP exit
CNVG? 00
SF .02

instead of the:

Code:

CNVG? 00
JMP exit

That's the method that is used in Lambert's W function as well as (sort of) in the Integrate routine. But I think we can save this additional iteration.

I set up a test program for the emulator and tested a dozen of different input values from 0,4999999999999999 to 1E-300, with 1 to 10000 degrees of freedom (including every single value between 1 and 100). The results showed an interesting pattern:

1. The number of uncertain (last) digits in a perfectly converged result is quite uniformly the number of digits in the degrees of freedom. In other words: for n=1...9 you can expect 33 valid digits, for 10...99 it's 32 digits, for n=100..999 it's 31 digits, and so on. I assume this is because of the term x² + n within the iteration loop.

2. This level of accuracy is reached as soon as the first 12 digits agree. So with 14 matching digits we should be on the safe side, as already mentioned. I think that this should be sufficient even for a 39-digit C function.

3. and most important: further iterations do not seem to guarantee a better result. Unless CDF(t) – p is exactly zero, they just cycle through various permutations of the last digits.

Example: p=0,1 with 20 degrees of freedom.

Initial guess for t: 1,312400768070987055677223263645284
After 1st iteration: 1,325340121782327835195769644183172
After 2nd iteration: 1,325340706985046343045588015602724
After 3rd iteration: 1,325340706985046343099953113344022
*** at this point CNVG 00 considers the quantile converged ***
After 4th iteration: 1,325340706985046343099953113344020
After 5th iteration: 1,325340706985046343099953113344022
After 6th iteration: 1,325340706985046343099953113344020
...

True 34-digit value: 1,325340706985046343099953113344018

4. There are cases where an additional iteration corrects the last digit(s) and then the approximation does not change any more. This does not neccessarily mean that the added iteration improves the result.

Example: p=0,05 with 20 degrees of freedom

After 3rd iteration: 1,724718242920787272832061725382084

This is the exact 34-digit result. Further iterations change this to ...082 and then stay there.

That's why I think that CNVG 00 should be fine here.

Dieter