[WP34s] Regularized incomplete Beta function

05082014, 11:29 PM
Post: #41




RE: [WP34s] Regularized incomplete Beta function
(05082014 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 

05092014, 09:06 PM
(This post was last modified: 05092014 09:09 PM by Dieter.)
Post: #42




RE: [WP34s] Regularized incomplete Beta function
(05082014 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. 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 62. 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 

05102014, 12:03 AM
Post: #43




RE: [WP34s] Regularized incomplete Beta function
(05092014 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 

05112014, 01:23 PM
Post: #44




RE: [WP34s] Regularized incomplete Beta function
(05102014 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 1E37. 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 1E6100. 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. 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 

05112014, 01:50 PM
(This post was last modified: 05112014 02:15 PM by Dieter.)
Post: #45




RE: [WP34s] Regularized incomplete Beta function
(05052014 08:26 AM)Paul Dale Wrote: After a bit of debugging, the new StudentT 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:
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 Tregister 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:
Dieter 

05112014, 09:51 PM
Post: #46




RE: [WP34s] Regularized incomplete Beta function
(05112014 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: The eight level stack is in play in xrom code. I just didn't document the stack diagrams fully.  Pauli 

05112014, 10:33 PM
Post: #47




RE: [WP34s] Regularized incomplete Beta function
(05112014 01:23 PM)Dieter Wrote: The CF branch uses a modified Lentz algorithm where the value for the "tiny" element is set to 1E37. 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 1E6100. What do you think? I've dropped the threshold to 1e10000. 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 

05122014, 06:05 AM
Post: #48




RE: [WP34s] Regularized incomplete Beta function
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 

05122014, 06:18 AM
Post: #49




RE: [WP34s] Regularized incomplete Beta function
(05112014 09:51 PM)Paul Dale Wrote:(05112014 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. 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 

05122014, 06:58 AM
Post: #50




RE: [WP34s] Regularized incomplete Beta function
(05122014 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 

05122014, 07:47 PM
(This post was last modified: 05122014 07:57 PM by Dieter.)
Post: #51




RE: [WP34s] Regularized incomplete Beta function
(05122014 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 Treplication or a simple rollup 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 

05122014, 08:14 PM
Post: #52




RE: [WP34s] Regularized incomplete Beta function
(05122014 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 Treplication or a simple rollup 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:) 

05122014, 09:59 PM
Post: #53




RE: [WP34s] Regularized incomplete Beta function
(05122014 07:47 PM)Dieter Wrote: Did you also say that the stack automatically is set to eight levels? This would make things like Treplication or a simple rollup 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 precomputed 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 

05222014, 08:51 PM
Post: #54




RE: [WP34s] Regularized incomplete Beta function
(05122014 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 

05252014, 09:43 PM
Post: #55




RE: [WP34s] Regularized incomplete Beta function
(05222014 08:51 PM)Dieter Wrote: In the meantime I have set up a similar routine for the Chi² quantile. Done. The results can now be found here. Dieter 

« Next Oldest  Next Newest »

User(s) browsing this thread: