[WP34s] Regularized incomplete Beta function

[Paul Dale]

(05-08-2014 10:58 PM)Dieter Wrote:  The user may enter dof values that lead to problems and strange results. At the moment I am looking at the Chi² quantile. For Chi² = n the CDF should be near 0,5 (slightly larger). Try this with n=2000 and you get 0,5042... (correct). Then try n=3000 and the 34s returns zero (!).

This is a bug, not a gradual loss of accuracy

The series expansion for the incomplete gamma function (as opposed to the continued fraction variation) isn't converging in 500 iterations for x around the degrees of freedom. For the n = 2000 case it is taking 440 iterations to converge. For n = 3000, it wants 532 iterations. For n = 250000, 4591 iterations are required.

Perhaps I got the threshold for switching between the continued fraction and the series expansion wrong.

This might also explain the gradual loss of digits you are seeing at time. That would match the behaviour of the continued fraction path which attempts to return something based on the best guess. The series code just returns zero when its iteration count expires.


Pauli

[Dieter]

(05-08-2014 11:29 PM)Paul Dale Wrote:  The series expansion for the incomplete gamma function (as opposed to the continued fraction variation) isn't converging in 500 iterations for x around the degrees of freedom.
...
Perhaps I got the threshold for switching between the continued fraction and the series expansion wrong.

I am not very good at finding the code of a particular 34s function on Sourceforge, so this is what I assume: the incomplete Gamma function may use the two "classic" series resp. CF representations as suggested e.g. in Numerical Recipies, chapter 6-2. Here, the recommended threshold for IΓ(a, x) is the point where x equals 1+a. For lower x the series is used, for higher x it's the continued fraction method. This would also explain why the Normal CDF requires significantly more time for x>√3, since the CDF is evaluated as IΓ(0.5, x²/2) so that the switch between both methods occurs here: 1 + 0.5 = (√3)²/2. If this really is the way IΓ works, you got the threshold right.

Quote:This might also explain the gradual loss of digits you are seeing at time. That would match the behaviour of the continued fraction path which attempts to return something based on the best guess. The series code just returns zero when its iteration count expires.

So does the code in the NR book. ;-)

I assume the CF is evaluated by some kind of Lentz's method. I do not know how fast this is in C, but in user code I prefer a different approach (not only) for speed reasons: estimate the number of required terms and evaluate the CF from right to left. I did this with the Normal CDF in a user code program for the 34s and it works very well and, most important, it's fast. I even think that this provides somewhat better accuracy. ;-)

Dieter

[Paul Dale]

(05-09-2014 09:06 PM)Dieter Wrote:  I am not very good at finding the code of a particular 34s function on Sourceforge, so this is what I assume:...

The function is decNumberGammap in decn.c line 2164. You assumptions are correct.

I've fixed the zero return to give a best guess return instead and increased the iteration counts which will help but won't avoid this issue.


Now to try to find a more robust algorithm. The GSL version has a third path for special argument values (for performance I suspect) but seems to use the same algorithm for the two main paths. It has an explicit failure case when it estimates that too many iterations will be required. I've seen the series and CF approaches used by other libraries too.


- Pauli

[Dieter]

(05-10-2014 12:03 AM)Paul Dale Wrote:  The function is decNumberGammap in decn.c line 2164.

Thank you. In my browser the line numbers on the left do not align with the code, but I see where the function code and the series resp. CF routines start.

The CF branch uses a modified Lentz algorithm where the value for the "tiny" element is set to 1E-37. This is a common value for single precision where the working range ends right below this. The examples in the NR book use a similar value since they obviously assume single precision as well. However, in our case I feel this value should be more like 1E-6100. What do you think?

Quote:I've fixed the zero return to give a best guess return instead and increased the iteration counts which will help but won't avoid this issue.
...
Now to try to find a more robust algorithm.

I do not have a solution yet. On the other hand I found some improvements for the Chi² quantile. The new code is running quite fast and it even handles cases where a simple Newton solver would fail. But it still needs some adjustments for high dof values with probabilities close to zero.

Dieter

[Dieter]

(05-05-2014 08:26 AM)Paul Dale Wrote:  After a bit of debugging, the new Student-T quantile function code appears to be working on both the 34S and the 31S. The next builds will get this

I see it is now in t.wp34s. Just two remarks:

1. Catching the case p=1/2 at the beginning is not neccessary. The estimate will be exactly zero, and this is what will be returned directly after the first iteration, i.e. almost immediately.

2. This is part of the iteration loop:

Code:


qf_q_Ix::   DROP
            x[^2]        // x^2 x x x
            ENTER[^]
            RCL+ J        // n+x^2 x^2 x x
            /
            Num 1/2
            RCL[times] J    // n/2 x^2/(n+x^2) x x
            Num 1/2
            [<->] ZXYX    // x^2/(n+x^2) 1/2 n/2 1/2
            I[sub-x]    // Ix 1/2 x x
            [times]
            +/-
            Num 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)

The stack shuffling command obviously was added because the incomplete Beta function now accepts the three parameters in a different order than before. After this, right before the call to the Beta function, the T-register holds 1/2. After that call, the comment says that both Z and T hold the current estimate x – which was completely lost before the call (was neither in X nor Y nor Z nor T). ?!?

I would say this should work better:

Code:


qf_q_Ix::   DROP
            Num 1/2        // 1/2   x    x    x
            ENTER[^]
            RCL[times] J   // n/2  1/2   x    x
            RCL T
            x[^2]          // x²   n/2  1/2   x
            RCL+ J         // x²+n n/2  1/2   x    (and x² in L)
            STO/ L
            X[<->] L       // x²/(x²+n)  n/2  1/2   x
            I[sub-x]       // Ix  x  x  x
            Num 1/2
            [times]
            +/-
            Num 1/2
            RCL- .00
            +

Dieter

[Paul Dale]

(05-11-2014 01:50 PM)Dieter Wrote:  1. Catching the case p=1/2 at the beginning is not neccessary. The estimate will be exactly zero, and this is what will be returned directly after the first iteration, i.e. almost immediately.

I'm getting "no root found" without that check. That's why I added it.

Quote:2. This is part of the iteration loop:

Code:

...

The stack shuffling command obviously was added because the incomplete Beta function now accepts the three parameters in a different order than before. After this, right before the call to the Beta function, the T-register holds 1/2. After that call, the comment says that both Z and T hold the current estimate x – which was completely lost before the call (was neither in X nor Y nor Z nor T). ?!?

The eight level stack is in play in xrom code. I just didn't document the stack diagrams fully.


- Pauli

[Paul Dale]

(05-11-2014 01:23 PM)Dieter Wrote:  The CF branch uses a modified Lentz algorithm where the value for the "tiny" element is set to 1E-37. This is a common value for single precision where the working range ends right below this. The examples in the NR book use a similar value since they obviously assume single precision as well. However, in our case I feel this value should be more like 1E-6100. What do you think?

I've dropped the threshold to 1e-10000. We're using the high precision internal reals here and exponents are very wide. I also noticed a bug where I didn't take the reciprocal properly. The Chi² problems you found for high degrees of freedom, aren't impacted by this change. The series approximation is used there.

Quote:I do not have a solution yet. On the other hand I found some improvements for the Chi² quantile. The new code is running quite fast and it even handles cases where a simple Newton solver would fail. But it still needs some adjustments for high dof values with probabilities close to zero.

Sounds useful of course.


Pauli

[Paul Dale]

Found a possibility in the GSL sources. For large n and x close to n, they switch to the continued fraction approach instead of the series. I've implemented something similar.


- Pauli

[Dieter]

(05-11-2014 09:51 PM)Paul Dale Wrote:  

(05-11-2014 01:50 PM)Dieter Wrote:  1. Catching the case p=1/2 at the beginning is not neccessary. The estimate will be exactly zero, and this is what will be returned directly after the first iteration, i.e. almost immediately.

I'm getting "no root found" without that check. That's why I added it.

Then something is going wrong. The first estimate should be zero, so the CDF is exactly 0.5 and the following correction zero as well. CNVG? 00 then compares two zeros in X and Y which should be considered converged. It works fine here with my original code on a hardware 34s. Please check.

BTW, if "no root found" appears, the value in X now is |quantile|, without a sign adjustment. That's why I originally used flag 01 to put the error message behind the final +/–.

Dieter

[Paul Dale]

(05-12-2014 06:18 AM)Dieter Wrote:  BTW, if "no root found" appears, the value in X now is |quantile|, without a sign adjustment. That's why I originally used flag 01 to put the error message behind the final +/–.

In XROM, the error causes the stack to revert to what it was before the operation commenced.

I don't understand why there is the problem with 1/2 but will look into when I get a chance.


- Pauli

[Dieter]

(05-12-2014 06:58 AM)Paul Dale Wrote:  In XROM, the error causes the stack to revert to what it was before the operation commenced.

I see. Probably that's due to the exit code that restores the stack after the XROM routine has finished. Obviously there are several things to consider when a regular user program is transferred to XROM. Did you also say that the stack automatically is set to eight levels? This would make things like T-replication or a simple roll-up command work differently. Is there some kind of list with all the details that are different in XROM?

My idea originally was to exit with an error message, but with the last approximation in X so that the user can decide whether to use this value or not.

Quote:I don't understand why there is the problem with 1/2 but will look into when I get a chance.

Just check the first guess for p=1/2. Is it zero? Then check the updated approximation. For p=0.5 the qf_q_Ix branch is used. The Beta function should return zero.

Dieter

[walter b]

(05-12-2014 07:47 PM)Dieter Wrote:  Obviously there are several things to consider when a regular user program is transferred to XROM. Did you also say that the stack automatically is set to eight levels? This would make things like T-replication or a simple roll-up command work differently. Is there some kind of list with all the details that are different in XROM?

Please see App. H of the manual.

d:-)

[Paul Dale]

(05-12-2014 07:47 PM)Dieter Wrote:  Did you also say that the stack automatically is set to eight levels? This would make things like T-replication or a simple roll-up command work differently. Is there some kind of list with all the details that are different in XROM?

On xIN, the stack goes to eight levels (but can be changed back to four if required), the device also enters double precision mode (which can't be changed back). These are the big two and xIN/xOUT is sufficiently useful to make putting up with them worthwhile.

The other changes are pretty minor. The A..D-> instruction becomes useful to create a small number of temporary registers but LocR is better practice. Branches lose their labels -- they are either replaced with relative branches or operate via a pre-computed jump table. CNVG? sometimes changes how it operates slightly to provide better single precision results.


The biggest change that causes most problems moving code to XROM is the eight level stack.

- Pauli

[Dieter]

(05-12-2014 06:58 AM)Paul Dale Wrote:  I don't understand why there is the problem with 1/2 but will look into when I get a chance.

Fine. In the meantime I have set up a similar routine for the Chi² quantile. I tried the same Halley approach as in the Student case, and it worked very well – until the number of dof got beyond 100 or 200. At this point either an extremely good estimate is required or the approximation will converge very slowly at first (until the approximation is very close to the true value and rapid convergence sets in). That's why I tried a different approach. The Normal quantile guess also got an update and now is more accurate than before. Only three constants need to be adjusted.

I think I'll post the results later in a new thread.

Dieter

[Dieter]

(05-22-2014 08:51 PM)Dieter Wrote:  In the meantime I have set up a similar routine for the Chi² quantile.
(...)
I think I'll post the results later in a new thread.

Done. The results can now be found here.

Dieter