[WP34s] A set of new quantile functions

05252014, 09:40 PM
(This post was last modified: 05282014 01:23 PM by Dieter.)
Post: #1




[WP34s] A set of new quantile functions
Some 34s quantile functions for the statistical distributions do not seem to work as fast and reliably as they could. I recently suggested an update for Student's t quantile, based on a new estimate and the Halley method for fast convergence.
I now would like to propose a similar solution for the Chi² quantile. Here things are a bit more complex: a regular Halley iteration can be set up easily, as shown in another earlier thread, but there are issues with large degrees of freedom. Although the first Chi² estimate is only approx. 11% low in the worst case, this can mean that the CDF of this can be off by several orders of magnitude. This leads to very slow convergence in the first 2...10 iterations which essentially apply the same correction until the estimate is sufficiently accurate and rapid convergence starts. So I chose a different approach: instead of cdf(x) – p, better solve ln(cdf(x)/p) = 0. This requires a bit of calculus, but it is not too hard to do. The Student case was already discussed earlier, so the code below is just a slight update. Please note that the incomplete Beta function recently was updated and now requires a new order of its three arguments in X, Y and Z. The code below assumes the previous order. Compared to the current version included in the 34s, the Normal quantile was streamlined and shortened. It also takes advantage of the new Normal estimate. The result below is part of a set of programs for the Normal, Student and Chi² quantile. They all rely on a short routine that provides an estimate for the Normal quantile. Compared to the routine previously used in the 34s, only some constants were changed to provide better accuracy. More important, "GNQ" now accepts any probability between 0 and 1 (exclusively) and it returns a correctly signed result. As an additional benefit, min(p, 1–p) is returned in Y which can be handy here and there. Finally, the following set of programs is what I got until now. This code is to be run in double precision mode, four stack levels are sufficient. As usual, the number of degrees of freedom is assumed in register J. Otherwise the programs only use R00 and R01 and flags 00 and 01. I could not find any severe errors, but please do your own thorough tests and do not consider this perfect. As usual, any error reports and suggestions are welcome. Code:
Edit: there was an update in the incomplete Beta function. If your 34s displays this function as "Iβ" you are running the older version and the original code above works fine. If the function is named I_{x} instead, this is the newer version with a different order of the three arguments. In this case use the following code between LBL 02 and 03 in "TQF": Code: ... And finally: Pauli, what do you think ?) Dieter 

05262014, 07:40 AM
(This post was last modified: 05262014 07:42 AM by walter b.)
Post: #2




RE: [WP34s] A set of new quantile functions
AFAICS, Pauli uploaded some modifications following your suggestions last night. Will be in next build of WP 34S and WP 31S.
d:) 

05262014, 08:06 AM
Post: #3




RE: [WP34s] A set of new quantile functions
These new distribution functions aren't in the image. I'll get to them once things calm down a bit here.
The changes were to the L.R. command.  Pauli 

05262014, 10:45 AM
Post: #4




RE: [WP34s] A set of new quantile functions
A new build will depend on Marcus. The L.R. changes are committed.
The distribution changes aren't going to happen immediately. I've more important things to worry about than calculators at present.  Pauli 

05262014, 06:43 PM
Post: #5




RE: [WP34s] A set of new quantile functions
(05262014 09:36 AM)fhub Wrote: That's a bit strange in my opinion: Yes. The essential point was not the order of the Iβ arguments, but the fact that the function worked different from what the manual said. This has been corrected. Iβ now works "as advertized", which indeed is more logical. But there is another point. The beauty of the German language has gifted us with the term "Die normative Kraft des Faktischen". What do you think, how many 34s are currently using the latest firmware? Five percent? Ten? That's why I assumed the previous Iβ implementation. All others may add a simple stack shuffle command. Quote:It would really be nice to have a (quite) bugfree WP34s version finally That's why I humbly proposed these three routines. The previous versions sometime threw errors or required inacceptable execution times. What about you? Your TVM solver has made it into the 34s and added some essential functionality. Great. But what about another project? While I worked on the quantile functions I discovered both the Iβ bug as well as the L.R. problem. Both have been corrected – and the 34s became a bit better. This is a community effort – "ask not what the 34s can do for you – ask what you can do for the 34s project". :) Quote:(and I don't care whether any highlevel function is accurate to 37 or 'only' 36 digits internally ). It's more like 30...34. There are limitations in the algorithm and the functions it uses as well as in the hardware. But I think 30 or 32 digits are good enough. If you want all 34 in DP mode, having these functions coded in C is the only way. ;) Dieter 

05262014, 06:52 PM
Post: #6




RE: [WP34s] A set of new quantile functions
(05262014 10:45 AM)Paul Dale Wrote: The distribution changes aren't going to happen immediately. I've more important things to worry about than calculators at present. Take your time. And most important: the suggested code is just a humble proposal. I do not say it is perfect and I strongly suggest you (and the rest of the community) do some thorough tests before this is included. Maybe there is a better approach, a faster solver, a better χ² estimate than good old WilsonHilferty, or whatever. I only can say that up to now everything seems to work fine here on my hardware 34s – within the limiations of the respective CDFs and Iβ resp. IΓ. Maybe these could get updated as well...?) Dieter 

05282014, 01:29 PM
Post: #7




RE: [WP34s] A set of new quantile functions
(05252014 09:40 PM)Dieter Wrote: Finally, the following set of programs is what I got until now. Addendum: I updated the original post and added a modified piece of code for use with the newest 34s firmware with the modified incomplete Beta function. I also noticed that this function now appears as I_{x}, which I think is a good idea. Dieter 

05282014, 02:19 PM
Post: #8




RE: [WP34s] A set of new quantile functions  
05292014, 12:43 PM
(This post was last modified: 05292014 02:12 PM by Dieter.)
Post: #9




RE: [WP34s] A set of new quantile functions
(05262014 10:45 AM)Paul Dale Wrote: The distribution changes aren't going to happen immediately. That's fine since another update may be required as well. I did some further tests with the recently suggested Student quantile function and found that close to zero the CDF of the calculated quantile often significantly differed from the original probability. Some tests with 42digit precision showed that the quantile was fine (33 or 34 correct digits), while there are some issues with the WP34s t(x) function. Near t=0 the Student CDF should be roughly 0,5 + r·t where r is a factor between 1/pi and 1/sqrt(2·pi), i.e. somewhere between 0,3 and 0,4. So the double precision results for the Student CDF should look like this: Code: t CDF(t) But actually there is a much larger t below which a plain 0,5 is returned: Code: t WP34s t(x) function for 10 dof For further investigations I used a little program that provides a 34digit tvalue of the same magnitude as entered (pi/3 * input): Code: LBL A Since the true result is approx. 0,5 the returned error can be read as "units in the xth digit". So for t near 5 E–17 the error is 2 units in the 17th digit. In other words, accuracy degrades down to single precision level. I assume this happens due to digit cancellation when the argument x = n/(t²+n) for the incomplete Beta function is evaluated, since for very small t this rounds to 1. That's why I suggest that for small t the Student CDF should be evaluated using 1–x = t²/(t²+n), as shown in the routine above. BTW the recently posted Student quantile function is not affected by all this as it already uses the suggested method for t < 0,5. EDIT: I looked at the current Student CDF code in t.wp34s, and indeed it works as assumed. Here is a fix that uses a different method for small t less than 1: Code: ... Dieter 

05302014, 06:13 AM
Post: #10




RE: [WP34s] A set of new quantile functions
(05252014 09:40 PM)Dieter Wrote: I now would like to propose a similar solution for the Chi² quantile. Here things are a bit more complex: a regular Halley iteration can be set up easily, as shown in another earlier thread, but there are issues with large degrees of freedom. Although the first Chi² estimate is only approx. 11% low in the worst case, this can mean that the CDF of this can be off by several orders of magnitude. This leads to very slow convergence in the first 2...10 iterations which essentially apply the same correction until the estimate is sufficiently accurate and rapid convergence starts. Recent tests showed that a similar behaviour may also occur with the Student quantile function. Up to now I only noticed this in extreme cases out of the SP working range, and here six or seven iterations seem to be sufficient. So the number of allowed iterations initially stored in R01 should be set to 8 for the moment. I will take a closer look at this and see if an approach similar to the Chi² case will eliminate this. Dieter 

06012014, 01:18 PM
Post: #11




RE: [WP34s] A set of new quantile functions
(05302014 06:13 AM)Dieter Wrote: Recent tests showed that a similar behaviour may also occur with the Student quantile function. It does. Here is an improved version of the Student quantile function that uses the same approach as the Chi² quantile: it solves ln(cdf/p) = 0. It does not require more than four stack levels. With eight levels the whole thing might get shortened a bit. The 34s code is rather complex and hard to understand, so here it is in a more readable form: Code:
The 34s code is the same as before, only the iteration loop at LBL 01 is replaced with this: Code: ... As before, 5 iterations are initially stored in R01. If the new code is too long and/or uses too much memory, the iteration limit in R01 may be removed. More than 5 iterations are only required in a limited domain for very small p and high dof. Up to now I have not found any cases within the SP working range. So if some more iterations are OK in these cases, the original code may be used as well. The worst I've seen so far are 12 loops (at 5000 dof, p=10^5400). However, the new code does it in three. ;) Dieter 

06022014, 08:32 PM
Post: #12




RE: [WP34s] A set of new quantile functions
It's me again. Here is the third reply to my own post. #)
(06012014 01:18 PM)Dieter Wrote:(05302014 06:13 AM)Dieter Wrote: Recent tests showed that a similar behaviour may also occur with the Student quantile function. Final update: there is another solution. The original (shorter) TQF version can be used with a slight modification. The program uses two different estimates. One for the central region and another for p close to 0 or 1. At least one of these seems to be good enough to ensure reasonably fast convergence. But in some cases, especially for very high dof values, choosing the one or the other is quite tricky since there is just a narrow corridor of probabilities where both approximations work equally well. So the crucial point is the decision whether to use the central or the tail estimate. Up to now this was done by checking whether the input is greater or less than 12^{–n}. This mostly works well, but it is not sufficiently exact for a few critical cases. After some tests I can suggest a different threshold that seems to work fine: simply replace the constant 12 with the term sqrt(n) + 7,5. So these lines near the beginning of TQF... Code: ... ...are replaced with these: Code: ... And that's it. So there are essentially two ways of solving the problem of slow convergence in some extreme cases: 1. The safe method: use a robust algorithm. This is the approach of the second version I posted, i.e. solving ln(cdf/p) = 0. However, this requires more code and memory. 2. Since the initial estimates are good enough, the threshold for switching between the two may be finetuned. With the suggested fix the original shorter version of TQF should work fine. I could not find a case where more than four iterations were required. Dieter 

06022014, 10:06 PM
Post: #13




RE: [WP34s] A set of new quantile functions  
06032014, 08:35 PM
(This post was last modified: 06032014 08:35 PM by Dieter.)
Post: #14




RE: [WP34s] A set of new quantile functions
(06022014 10:06 PM)Paul Dale Wrote: I'm reading these, just not acting at present. Thank you for your feedback. Today I found an even better way of defining the threshold that works far more accurately than required. Should this be used anyway? Yes, I think so. Why? Because we can. ;) So if the updated Student quantile should eventually make it into the 34s, please let me know. But again: take your time. Dieter 

« Next Oldest  Next Newest »

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