![]() |
[WP34S] Curious Bug in Inverse Normal Function - Printable Version +- HP Forums (https://archived.hpcalc.org/museumforum) +-- Forum: HP Museum Forums (https://archived.hpcalc.org/museumforum/forum-1.html) +--- Forum: Old HP Forum Archives (https://archived.hpcalc.org/museumforum/forum-2.html) +--- Thread: [WP34S] Curious Bug in Inverse Normal Function (/thread-219710.html) |
[WP34S] Curious Bug in Inverse Normal Function - Les Wright - 05-01-2012 The XROM based function is perceptibly slower, which is a fair trade given the huge gains made in extra flash memory by getting this and other functions out of the native code. However, in testing numerous arguments, I stumbled across one that creates a hang, for no good reason--namely, exactly p=0.0025, as well as its complement p=0.9975. Running the function on either of these probabilities makes the calc hang, requiring a reset. However, adding even the smallest perturbation--say 1e-16--allows function to converge without difficulty. Weird. Any idea why these particular values prove to be so disagreeable? Also, for what it's worth, I have written my own inverse normal routine that is at least as fast if not faster than the present built-in function, and which handles this quirky probability without problem. I wonder if the inverse normal code got broken or at least bruised in its move to XROM?
Les
Re: [WP34S] Curious Bug in Inverse Normal Function - Paul Dale - 05-01-2012 No idea what is wrong here but I can definitely reproduce it. Will look into it later. I'd be interested in your inverse normal implementation if it can handle the full range of values and gives a sufficient number of digits of accurate result. e.g. what does it do given an argument of 1e-350 ? .999999999999 ? 0.5 ?
Re: [WP34S] Curious Bug in Inverse Normal Function - Paul Dale - 05-01-2012 I know the problem. The iterative process to refine the result is oscillating between two values:
-2.807033768343804117221810394712332 The old inverse normal distribution code carried extra digits and this wouldn't have been an issue. Now to figure out a fix.
Re: [WP34S] Curious Bug in Inverse Normal Function - Paul Dale - 05-01-2012 Try now :)
- Pauli
Re: [WP34S] Curious Bug in Inverse Normal Function - Les Wright - 05-01-2012 Quote: No and yes. I haven't included a test that reflects the sign of the result for arguments above 0.5, but that is no problem for my purposes. It handles the first and third examples you give very well. The second I would reflect to 1e-12 and switch the sign. It is basically a stripped down version of your own routine--I get the initial estimate using the formulae Dieter provided, and simply do three successive Halley refinements against the calc's Phi(x) whether I need them or not. No convergence tests. As you know, two Halley refinements give at least 18 digits (in theory), and three give at least 54 in an environment of adequate intermediate precision. More than enough for the calc's SP and DP. For SP it is overkill, but it is still fast enough for me. I suspect I lose a digit here and there due to rounding, as I can't take advantage of the full 39 internal digits, but it really was a fruitful programming exercise for me. Here is my listing--forgive me for not commenting on it due to being rushed:
0001 ****LBL'INN' My sense is that you would do well in your XROM version just to do three Halley refinements and not bother with tests of convergence. It may seem overkill, but I bet it will speed things up.
Re: [WP34S] Curious Bug in Inverse Normal Function - Les Wright - 05-01-2012 Quote: The last five digits of the true result, properly rounded, would be 12269. Your fixed version goes with the former of the above two choices, 12332. My code, as listed above, gets a little closer when run in DP, with the last five being 12308. More than close enough for real life :)
Les
Re: [WP34S] Curious Bug in Inverse Normal Function - Marcus von Cube, Germany - 05-01-2012 A probable solution (instead of just relaxing the convergence criteria) could be to check for a cycle (keep two or three guesses and check for duplicates) and if one is found, disturb the value by one or a few ULPs, and try again. Do any iteration only a finite number of times.
Re: [WP34S] Curious Bug in Inverse Normal Function - Dieter - 05-01-2012 Ah, the Normal quantile again. :-) I would like to add some remarks that might help a bit here.
The applied algorithm itself is quite stable - as long as there is sufficient precision both in the CDF and, even more important, while evaluating the term CDF(x) - p This expression is used both in the Halley method as well in another approach I tend to prefer. The more x approaches the true quantile, i.e. the better CDF(x) matches p, the more important this point gets as the number of significant digits in the difference decreases. Another problem is the fact that, at least for x < 1, there are multiple n-digit values for x that will return the same n-digit CDF, so the nominator may become zero even if x is not yet exact.
My current solution (used in a 41 and 35s) works as follows:
i.e. p = 0,1 => x = +1,28155... Initial estimate, using a very simple rational approximation:The refinement is based on the method suggested in Abramovitz & Stegun, p. 954. I think it works very well. ;-) I originally used this method for the 41C with its 10-digit precision and a range down to 1E-99. With at least 12 digits working accuracy, one single iteration step should provide 11 valid digits +/- 1 ULP, i.e. it should be fine for a machine code routine that carries 13 internal digits before it returns an exact 10-digit result to the user. During my tests in Excel (15-digit binary arithmetics) I did not notice any error beyond +0 and -4 units in the 12th significant digit.
For calculators with a working range down to 1E-500 a different rational approximation for the initial estimate may be used, or the latter gets a slight adjustment: if x > 21,3 then x = x - exp(x/21,3 - 7,7)Alternatively, one or two more iterations are fine, too - and very fast as well since the PDF is evaluated very quickly far out in the distribution tails. The desired 16-digit accuracy for the 34s should be no problem with just two iterations, or three in double precision mode. That's why I agree that a fixed number of iterations may avoid the problem of oscillating approximations while the desired accuracy is still achieved. On the other hand I would recommend a different way of calculating the quotient t. For instance, the ratio CDF(x)/PDF(x) usually shows up during the calculation of the CDF (by series or continued fraction), so this might be used for an alternative way to calculate t.
Dieter
Re: [WP34S] Curious Bug in Inverse Normal Function - Dieter - 05-01-2012 Les, I think it's simply impossible to get a result with n exact digits as long as the working precision here is just the same n digits. ;-) Try it with a simple 10-digit machine - if you get 9 exact digits +/- 1 ULP you've achieved the best that can be done. That's why there are guard digits. ;-)
So if the working precision is limited to 34 digits, a result with 30-32 valid digits is fine. And that's what we get here - simply round the result to 32 digits and we're done: x = -2.8070 33768 34380 41172 21810 39471 23Dieter Re: [WP34S] Curious Bug in Inverse Normal Function - Paul Dale - 05-01-2012 I'm thinking using the same solver as the other distributions might work better here. It handles these cycles cases (switching to a bisection instead) but is going to be a bit slower.
- Pauli
Re: [WP34S] Curious Bug in Inverse Normal Function - Paul Dale - 05-01-2012 Quote: This is fundamentally the problem. Moving these functions to xrom means I've lost five digits of working precision, their accuracy must suffer because of this. It clearly suffered more than I'd thought :-(
Re: [WP34S] Curious Bug in Inverse Normal Function - Dieter - 05-01-2012 Just for the record: I did a quick-and-dirty implementation of the algorithm I proposed earlier in this thread, and I tried the given example p = 0,0025. Et voila - after just two iterations in double precision mode... Quote:...this is exactly what I got: x = 2,8070 33768 34380 41172 21810 39471 2269This does not change in further iterations. And yes, I know, this example does not prove anything. #-)
Dieter
Re: [WP34S] Curious Bug in Inverse Normal Function - Marcus von Cube, Germany - 05-02-2012 Pauli will be listening, I'm sure! We should get the normal distribution right because it's a cornerstone function. If we do get the same level of accuracy out of the other distributions, fine, if not, just force a rounding to 30 digits (in double precision mode only, of course!).
Dieter, your knowledge about the numerical algorithms for the probability functions is outstanding! Where does it come from? Are you a mathematician or just a genius? ;-)
Re: [WP34S] Curious Bug in Inverse Normal Function - Paul Dale - 05-02-2012 I'm listening. If I get keen enough, I'll give this method a trial -- it will cost more bytes in flash but will be faster and likely at least as accurate. At the moment, I'm a bit burned out over coding statistical functions as keystroke programs.
-2.807033768343804117221810394712233 Slightly better but not perfect. I doubt we'll be able to get much better without carrying extra digits.
Re: [WP34S] Curious Bug in Inverse Normal Function - Paul Dale - 05-02-2012 I think I am going to have to revisit the inverse normal function yet again. Any chance you could send the quick and dirty code please? Even better a commented version :-)
- Pauli
Re: [WP34S] Curious Bug in Inverse Normal Function - Les Wright - 05-02-2012 Quote: In any case, I think it will be wise to heed Dieter's warnings and be especially attentive to the behaviour around p = 0.5, where most digit loss is bound to occur due to that phenomenon of more than one n-digit quantile estimate associated with the same n-digit probability. Indeed, the distribution-formulas readme accounts discusses this in one of the suggested approaches, and does something special when p is in the interval (.495, .505).
A savvy user who is really obsessed with high accuracy (like me, even though that serves little practical purpose :) ), will be aware of the relationship between the inverse normal and the inverse chi-square with one df. Basically, for p < .5: Phi(-1)(p) = -Sqrt(Chi^2INV(1df)(1-2p))
For example, for p = 0.4999, the Phi(-1)(x) built-in function gives: -2.506628300880350989206501054077314 e-4
In comparison, passing 0.0002 to Chi^2INV(x) with 1 in register J and taking the negative square root yields: -2.506628300880350989206501054078343 e-4 The direct method rounds to 30 good digits. Using the inverse chisquare intermediary yields ALL 34 digits correctly.
About the worst example I have been able to find is p = .4999999, where the built-in function gives: -2.506628274631026751765674822177285 e-7 The inverse chi-square path gives:
-2.506628274631026751765674822754539 e-7 Again, the latter gives full 34-digit accuracy. The former is only 28 digits to 1 ULP, or 29 digits to 6 ULP. Just food for thought.
Re: [WP34S] Curious Bug in Inverse Normal Function - Dieter - 05-03-2012 Pauli, Les et al. - I'll try to summarize all this in one single post. ;-) Yesterday I replaced the initial rational approximation from by previous suggestion by a dedicated version for p down to 1E-400. It now provides an estimate with an absolute error near +/- 0,002 and requires just two short real constants. The rest can be done with INC commands. ;-)
Estimate for 1E-400 =< p =< 0.5: r = sqrt(-2 ln p)An estimate x below 0.001 should get rounded to zero to ensure best results in the following correction step. The approximation was designed to give eight valid digits after the first iteration (max. error approx. +/- 2 units in the 9th digit), and with a second resp. third iteration we should be on the safe side for standard resp. double precision - as long as the working precision is sufficient. But that's another story - see my other remarks below. Here's my implementation for the 34s.
001 **LBL AThis works fine at least for quantiles >= 1. Now let's look at p close to 0,5. Note: "+1D34" = +1 unit in the 34th significant digit. 0,49 [A] -2,5068 90825 87110 35762 36343 18346 9043 E-2 +1D34This shows the problem Les and I have been talking about previously. Generally spoken, our basic approach (i.e. evaluate the CDF for the current approximation, see how good it matches p and apply a correction) is good for up to 34 decimals, but not for 34 significant digits. This is because close to zero there are numerous n-digit quantiles that return the desired p so that no further correction is applied. In other words, a calculated quantile near zero with a magnitude of 10^-n may be off in its last n digits, so that we may expect something like 34-n digit accuracy. Yes, the given examples are one digit better. ;-)
What can we do? We will have to handle p close to 0,5 separately. This can be accomplished in several ways:
But there's also another way to address this issue: Forget about 34-digit accuracy. There simply is no way to ensure a valid n-digit result as long as all calculations carry just the same n digits and no further guard digits. I remember how I tried to get an exact 10-digit result for the Normal quantile on the 41C with the same 10 digits working precision. It's simply not possible (at least I do not know of a method to accomplish this). That's why I designed the algorithm I suggested a few days ago - it returns 11 good digits over the 41C's working range as long as the calculation is done with a few more guard digits, e.g. in machine code. If I remember correctly, the original idea behind the double precision mode was this: enable the user to write user-code functions with lots of extended precision, so that the final result can be assumed to carry at least the usual 16 digits that are normally returned to the user. Once again: in many, if not most cases it's simply not possible to get a true 34-digit result unless the working precision exceeds these 34 digits. So we cannot expect 34 digits anyway. What is possible? 32? 30? 25? It depends. Round it to 16 digits and you are on the safe side. :-)
Dieter
Re: [WP34S] Curious Bug in Inverse Normal Function - Paul Dale - 05-03-2012 Quote: Thanks very much for this, I'll have a look at including it or something similar. We've not so constant adverse any more. Constants for C occupy 36 bytes plus the code to load them, for keystroke programs they require nine bytes for single precision and 17 for double (plus two bytes to load them). Additionally, integers 0 - 255 are directly available and with the decimal point shifts, many short constants are only two instructions (4 bytes) to get.
Re: [WP34S] Curious Bug in Inverse Normal Function - Dieter - 05-04-2012 Hmhmhmhmhmmmm.... ;-)
I tried your examples with the emulator. It's a somewhat older version (3.0 2294). Here is what I got with the internal function, i.e. based on 39 digit internal precision: p = 0,49So the internal function seems to work quite well (with the one exception at 0,49999 999). I also compared the results to those obtained with the Chi-square method, but the latter often showed larger errors. For example: p = 0,49999 999What's going on here? Why do we get different results? Is it because my version still uses the internal 39-digit precision and yours is based on 34-digit XROM code? Dieter
Edited: 4 May 2012, 3:10 p.m.
Re: [WP34S] Curious Bug in Inverse Normal Function - Paul Dale - 05-04-2012 Most likely the extra digits being carried makes for the difference. Five bonus digits is a lot. There is also a small change to the Newton search which might also make a difference. The older code did some range limiting on Newton steps within the known interval of the solution) and reverted to bisection if things didn't seem to be converging fast enough, whereas the newer code doesn't do this. Bisection now only happens if the Newton step would go outside the known interval of convergence. The idea was to try to trap cases when convergence via Newton steps is very slow because each step only moved a small amount.
Re: [WP34S] Curious Bug in Inverse Normal Function - Marcus von Cube, Germany - 05-05-2012 Are phi(x) and PHI(x) exact to 34 significant digits? If not, it might help to leave them in as C code to improve the accuracy. I've no plans to implement quad precision in XROM. ;-)
Re: [WP34S] Curious Bug in Inverse Normal Function - Paul Dale - 05-05-2012 phi(x) will be close. PHI(x) several digits off. PHI^-1(x) even less so. Pushing these back into C is always an option.
Re: [WP34S] Curious Bug in Inverse Normal Function - Dieter - 05-08-2012 Pauli, I've been thinking a bit about the possible accuracy of the Normal quantile function. The 39-digit C-coded internal function seems to work fine for most cases, but there are a few flaws here and there, especially close to p = 0.5 resp. z = 0. In double precision user code (34 digits working precision) the result may even be wrong in the last six (!) digits. But we can do better.
Dieter
Edited: 8 May 2012, 7:40 p.m.
Re: [WP34S] Curious Bug in Inverse Normal Function - Les Wright - 05-08-2012 Dieter, can you take a closer look? It seems that your results for "Corresponding CDF" actually give more 9s than that specified in the desired probability, so I think there are typos here. If there are not, then I don't follow your analysis...
Re: [WP34S] Curious Bug in Inverse Normal Function - Dieter - 05-08-2012 Yes, sorry, the returned CDFs originally had a few nines too much here and there. #-)
Dieter
Re: [WP34S] Curious Bug in Inverse Normal Function - Les Wright - 05-08-2012 Never mind! According to the time stamps, you fixed the errors just minutes before I posted... ;)
Re: [WP34S] Curious Bug in Inverse Normal Function - Les Wright - 05-08-2012 Dieter, this is excellent! I can't wait to find the time to code it myself. Just to clarify things for myself--it seems you are dealing actually with the upper tail probabilities and confining your code to positive quantiles and probabilities less than 0.5, so according to convention we are talking about q, not p. I think this is really important to observe in the final code. If you want to make yourself cringe, compute, in double precision, inverse Phi for unity less 1 ULP (i.e, 1 ENTER X.FCN->ULP -). The result is good only to 3 or 4 digits. However, if you compute the quantile associated with the complement 1e-33, the result is almost exact to 34 digits (but with opposite sign, of course). Basically, the ultimate qf code really needs to take the complement of probabilities above 0.5, use that for the routine, and flip the sign. I also really think that the refinement should go back to Halley's method here, since in this case it is just so simple. Using your variables, dx(t) = t/(1-xt/2). What could be easier?
Les Edited: 8 May 2012, 10:59 p.m.
Re: [WP34S] Curious Bug in Inverse Normal Function - Paul Dale - 05-08-2012 Quote: Care to code it for xrom inclusion?
Re: [WP34S] Curious Bug in Inverse Normal Function - Paul Dale - 05-08-2012 Quote: Or support both p and q for the normal and cross call between them as required.
Re: [WP34S] Curious Bug in Inverse Normal Function - Dieter - 05-09-2012 Les, all algorithms I proposed for the Normal CDF and its inverse use the same convention: x >= 0 and p (or q) < 0,5. So, actually only the right half of the function (x is non-negative) is considered, and the upper tail integral (from x to infinity) is evaluated. This means:
In other words: at least my own CDF and quantile routines already do what you suggested. So the example you gave (quantile for 1 - 1E-33) actually is evaluated as -quantile(1E-33). I also proposed this about a year ago when the distribution functions were discussed here. I assumed this went into the actual implementation, but to my astonishment this obviously is not the case. Pauli, did something change here? Coded correctly, quantile(p) and quantile(1-p) should only differ in their signs, while they otherwise agree in all digits. So, yes, this idea of "the ultimate QF code" has been proposed. I just wonder why it is not implemented, at least not in the current version.
Regarding Halley's method: yes, it converges quite fast, but there is an even better approach (cf. Abramowitz & Stegun): Halley: dx = t / (1 - x*t/2)For small x*t/2 this approaches dx ~ t * (1 + x*t/2)Now the method I proposed adds one more term to this: A & S: dx = t + t^2*x/2 + t^3*(2*x^2+1)/6This approach converges a bit faster, and there is another advantage as well: for x = 0 it is equivalent to the sum of the first terms of the Normal quantile's series expansion. I used this nice feature here and there to provide exact results close to x=0 with the initial guess set to zero. By the way, my current double precision version for the 34s starts with the good old Normal estimate we discussed some time ago. So the rational approximation I recently suggested is not used here. This is because of the extremely wide working range in double precision mode that goes even below 1E-6000. The rational approximation works fine down to 1E-400, but farther out in the tails the usual three refinements would not suffice. That's why I reverted to the previous method, which is approximately asymptotic as p approaches zero. I just added a slight improvement: the term 0,2/u resp. 1/5u was replaced by 1/4u. The switch between the center and tail approximation is between 0,232 and 0,233, limiting the maximum absolute error to slightly less than 0,01. :-) Now I am looking at the CDF algorithm for large x. Currently I am using the internal function. The continued fraction loop is evaluated best (and fastest) from the right to the left, so the number of loops n has to be known beforehand. Currently this is done by an estimate for n I proposed last year. However, this has not been tested for 34+ digit precision. I also thought about the Lentz method so that the continued fraction can be evaluated from left to right until convergence is reached, but I feel this is way too slow here. Dieter
Edited: 9 May 2012, 7:37 a.m.
Re: [WP34S] Curious Bug in Inverse Normal Function - Dieter - 05-09-2012 I thought this already got implemented a long time ago?! Did something change here? To ensure the best possible accuracy, both the CDF and the QF should work the way I suggested in my reply to Les. Please see below.
Dieter
Re: [WP34S] Curious Bug in Inverse Normal Function - Les Wright - 05-09-2012 I use the Lentz method in my own adaptation of the right-sided incomplete gamma, and I do agree that when I compute erfc or the upper-tail normal as a special case of incomplete gamma it seems slower than it needs to be. Also, because the arguments have to be manipulated before being passed to my incomplete gamma routine, that's another potential source of digit loss due to rounding. You aspire to full internal 39-digit accuracy of the CDF (or, conversely, the upper-tail probability), but this is impossible with only 39 digits since, even if every single computation returned 39 perfect digits and there was no loss due to rounding there is the matter of subtraction error at the threshold where one moves from the series expansion to the continued fraction. You have chosen 2.326 as your break point, and I think this is a good one, since that is where the upper tail probability is is 0.01, but even as you approach that you will lose at least one digit to subtraction of the series result from 0.5 to yield the upper-tailed probability. If you make the threshold lower, digit loss is reduced, but the continued fraction gets much more sluggish as more--sometimes many more--terms are needed and that means more computation and accumulated rounding errors. The existing code uses a threshold of roughly 3.15. The series does quite well up there, and beyond it the continued fraction is very fast. However, the subtraction error is at least two or three digits lost when we subtract the series result from 0.5 to get the upper tail probability. If we push the normal CDF and upper-tail probability and PDF back to C code, 34 digits returned to the stack should be easy to attain, but if we insist on 39 digits internally the number type used will have to introduce some more guard digits. That subtraction error is a bugger, but your use of the 2.326 threshold is probably the best compromise between speed and accuracy. Whether due to subtraction or rounding these calculations will not across the range give 39 digits of internal precision all the time unless working precision is more. However, the 34 digits we see on the stack in DP should be easy.
Les Re: [WP34S] Curious Bug in Inverse Normal Function - Dieter - 05-09-2012 Les, I know that it is not possible to get 39-digit accuracy on a 39-digit machine. This was already pointed out in several earlier posts, so on my 12-digit 35s I am happy with 11 valid digits (+/- 1 unit) in the final result. That's why, on the 34s, I would like to see the CDF (and, if possible, also the QF) implemented in an internal function with 39-digit working precision - to achieve a final result with 34-digit accuracy. ;-) Yes, switching from the series expansion to the CF method at x = 2,326 will sacrifice one digit. The ideal break point would be x = 1,28155... where Q(x) is 0,1. But this would slow down the calculation significantly. That lost digit is also the reason why I think the CDF (and quantile) should be coded in C with 39-digit precision (again: in order to get at least 34 valid digits in the result - we will not get all 39 anyway). I think we agree that the current break point at Pi or 3,15 is not perfect if the goal is maximum accuracy. This point is even beyond x = 3,09023 (Q(x) = 0,001) so we will lose even four digits, if I get it right. Maybe Pauli can change the code to switch at least x = 3 (will save one digit) or, even better, 2,326 (two more). On a side note, choosing 2,326 as the break point was the obvious option for my earlier 35s programs, since near this value the execution times for the series expansion and the continued fraction were roughly the same. ;-) And losing one digit was okay for me since a 12-digit machine will not produce a 12-digit result anyway. #-) It was and is exactly what you said: a compromise between speed and accuracy. In general, I think we both see the same problems as far as numeric accuracy is concerned. I suggested the one or other method to overcome these. If we assume that on an n-digit machine we may achieve about n-2 valid digits in the computed quantile (provided the CDF is exact!) this again leads to my initial suggestion: let's have these two key functions coded in C so that the result will have 34 safe digits. I just see that's exactly what you wrote in your last sentence. ;-) So we seem to agree in all major points on the numeric side. Now let's try to convince Pauli to move these two essential functions back to C. ;-)
Dieter
Re: [WP34S] Curious Bug in Inverse Normal Function - Paul Dale - 05-09-2012 Changing the break point in the xrom implementation is easy and cheap. Changing it in C is possible but requires more memory for the task. At the moment we've got the memory free to make some improvement to this function.
Quote: I'm happy to do some relocating. What functions exactly? The standard normal functions as they previously existed are 856 bytes in size -- minus the space they currently use. This is probably a bit on the large size when the current functions are returning 16 digit accurate values albeit rather slowly.
The QF could then be keystroke programmed.
Re: [WP34S] Curious Bug in Inverse Normal Function - Paul Dale - 05-09-2012 I can use more than 39 digits internally too -- some complex operations are performed carrying 72 digits, the modulo two pi carries a completely insane number of digits and the standard remainder operation uses higher precision due to a flaw in the decNumber library.
Re: [WP34S] Curious Bug in Inverse Normal Function - Paul Dale - 05-10-2012 I've moved the normal pdf back into C. It takes x in X and returns the pdf in X, 'a' in Y and 'b' in Z. Registers Y and Z are destroyed and last X is not set -- this is a very internal function. After executing cdfqh, you should execute * + to get the lower tail CDF result. Better is the sequence cdfqh 1 CPX_* DROP gets the same result but with an atomic multiply-add instruction that handles double precision results without loss of accuracy. That's most of our flash headroom gone now :-(
- Pauli
Edited: 10 May 2012, 12:49 a.m.
Re: [WP34S] Curious Bug in Inverse Normal Function - Marcus von Cube, Germany - 05-10-2012 Pauli, is this code conditionally compiled? Not everybody needs 34 digits for his statistical calculations.
Re: [WP34S] Curious Bug in Inverse Normal Function - Paul Dale - 05-10-2012 No it isn't a conditional include. This replacement code will be significantly faster than the existing xrom code and that justifies it place enough for me. The conversion to xrom slowed all the statistical code down a bit, the normal distribution suffered this by far the worst. More importantly, I don't want to have to support dual versions of the same thing :-) The quantile function won't be 34 digits accurate -- that simply isn't achievable from xrom.
Re: [WP34S] Curious Bug in Inverse Normal Function - Les Wright - 05-10-2012 For Phi(x) and Phi-sub-u(x), the constituent parts from the cdfqh helper function are returned to the stack and lose their guard digits, so there still is indeed the likelihood that things will be off a little often in the 34th digit--not the full 34-digit accuracy I know Dieter was hoping for. But, at least below the new threshold of |x| < 2.326, things seem fast and really good. But there seems to be a more serious challenge. You have implemented the new cutoff as noted. The last four digits of of the 34-digit Phi-upper-(2.326 - 1ULP) result are 7055, whereas the true result properly rounded is ends with 7050. Excellent 34-digit precision considering that the last few steps of the computation are done in 34-digit working precision, not 39. But when we get into continued fraction territory, there be wonkiness. The calculator's 34-digit output for Phi-upper-(2.326+1ULP) ends with 3078437, whereas the true last seven digits are 0377044. The may be evidence that the formula used to compute the number of terms taken with the CF needs to be reconsidered, since in comparing a calc with the newest firmware and one I haven't yet re-flashed it seems this is old news. For example, when I compute the upper-tail Phi(x) for 3.2, an argument well past the OLD threshold of Pi, the 34-digit output is "wrong" in the last 8 digits in both 2965 and 2975. Moreover, that "wrongness" agrees in all of those digits except the last. On the other hand, my slower and less efficient routine that computes the upper-tail probability for such larger arguments as a special case of the right-sided incomplete gamma function (my own routine) gets 32 of the 34 digits correct BOTH with the threshold examples and the argument of 3.2. I am sheepish I didn't pick up on this before, but I took for granted that the CDF and upper CDF were pretty solid in double precision whether in C-code or XROM. It may be worth revisiting whether a modified Lentz computation is worthwhile after all. Keep in mind the calculator's code already uses this algorithm for the betacf and gcf subroutines. Phi and Phi-upper are really just special cases of the left and right gamma function respectively, so maybe we should forgo a dedicated brute force right to left computation of the continued fraction and use the modified Lentz code that is there? I admit this would mean another byte-gobbling constant to pass to gcf() as the gln parameter (ln(sqrt(Pi)) = ln(Pi)/2 = 0.572364942924700087071713675676529355824), but if we ditch the dedicated CF code, would it be worth it? My user code of the right-sided incomplete gamma is perceptibly slower, but if I am getting 32 or 34 digits in a 34-digit environment, there is something right about this alternative. Just my 3 euros...
Les Edited: 10 May 2012, 3:12 a.m. after one or more responses were posted
Re: [WP34S] Curious Bug in Inverse Normal Function - Paul Dale - 05-10-2012 There are quite possibly errors in the code still :-( I started using gammap for the normal cdf but cancellation limited the number of correct digits at times -- this is what Dieter first picked up on. I really need to implement the right sides of the incomplete gamma and beta functions and use them for the other distributions (the upper tails for F, chi^2, binomial and Poisson are pathetic at present). User code for both sides of the incomplete beta and gamma would probably make things space neutral -- xrom seems to be about half the size of C for this kind of thing. Speed has also been subservient to accuracy during the development of the 34S.
PS: the upper tail CDF's were introduced this past week -- so you couldn't have taken them for granted.
Re: [WP34S] Curious Bug in Inverse Normal Function - Les Wright - 05-10-2012 Quote: Ah, but in a fashion it has always been there for the normal CDF. I have simply flipped the sign on the argument to get the upper CDF. Les
Edited: 10 May 2012, 4:24 a.m.
Re: [WP34S] Curious Bug in Inverse Normal Function - Dieter - 05-10-2012 Les wrote: Quote:I have not tried the new version yet (all I have is the emulator) but if it's not worse than my user code solution, we should get close to 34 digits. I think an error within 1 unit of the 33th digit is feasible. Quote:I tried these examples with an older emulator version (39-digit C) and in all cases the exact 34-digit value was returned: 2,326 [ENTER] ULP [-] => 0,010009 27534 08676 67260 76355 29037 7050That's why I prefer the CDF in C.
Then I tried the same values examples with my own 34s user code implementation, and this is what I got: 2,326 [ENTER] ULP [-] => 0,010009 27534 08676 67260 76355 29037 7049The true values are Phi(2,326 - 1E-33) = 0,010009 27534 08676 67260 76355 29037 70496 885...So in all cases the results are correct within 1 unit of the 34th digit. Quote:Yes, definitely. As I said before, the original estimate for the number of loops has been tested only for up to 16 digit accuracy. According to the original formula, the nominator in the estimate for the required number of loops would be something like 500 here, i.e. n = 500/(x-1) + a few to be safe. In fact it seems to me (just tried a few examples) that this number of loops is good for 34 digits, but with some error in the last place due to roundoff errors. We should be aware that the number of iterations must be as high as required, but it should be as low as possible to prevent excessive roundoff. So I tried n = 400/(x-1) + 8 (caveat - this has not been tested!) and I got the results you see above. I do not know how the number of iterations it determined currently, but I think I remember something like 200/(x-1), i.e. roughly half of what I tried. Which means that we will not get 30+ digits at all. And that's what you see in the results you posted. But I would even go one step further. I agree with Marcus in that not everyone needs 34-digit accuracy here. My suggestion is: let's make the special functions that require quite a bit of execution time accuracy-dependent. In other words: they should be evaluated to 16+ digits if this is what is returned to the user, or full 39 digits if we are in double precision mode, or that function is called internally by another one that in turn has to return a 34+ digit result. There is no need to wait for a 34-digit result if the user only gets 16 digits returned. The former would require about 4x the time of the latter. Compared to the current implementation this will roughly double (!) its execution speed, since for 16 digits around 100/(x-1) loops (instead of 200/(x-1)) are sufficent.
Dieter
Re: [WP34S] Curious Bug in Inverse Normal Function - Les Wright - 05-10-2012 Quote: Of course you would! In an older version you would be dealing with the series still for arguments around 2.325, since the cut-off was around Pi in that older code. I would expect the results would be pretty good, as you see. If you want to see how the continued fraction performed, you have to look at a larger argument. Try 3.15 or 3.2. I am sure you will see that there are not enough CF loops for 34 digits and you'll get 27 to 30 correct. I got a nice email from Pauli asking for an actual code sample, and I will see what I can do. In the meantime, I am starting to wonder if my obsession with full double precision accuracy for these functions is reasonable or practical. Getting even sixteen digits, only 12 of which are routinely displayed is pretty much all most reasonable humans need!
Les Edited: 10 May 2012, 3:31 p.m. after one or more responses were posted
Re: [WP34S] Curious Bug in Inverse Normal Function - Marcus von Cube, Germany - 05-10-2012
Quote:I doubt we are reasonable humans. ;-) Re: [WP34S] Curious Bug in Inverse Normal Function - Marcus von Cube, Germany - 05-10-2012
Quote:If I read Pauli's code correctly, the loop count is fixed at 500, not depending on X. It breaks out of the loop if the iteration does not change the value. Edit: I misread the code, the 500 loop is not for the continued fraction but for the series expansion below the threshold. The CF count is something like 256/(|x|-1). I might be able to make this dependent on the requested accuracy.
if (dn_lt(&absx, &const_2_326)) {
Edited: 10 May 2012, 12:58 p.m.
Re: [WP34S] Curious Bug in Inverse Normal Function - Dieter - 05-10-2012 Les, of course I expected the internal function to be exact. The point here was that the user-coded function also may return a result with 33 or maybe even 34 digits - provided the number of terms in the CF loop is sufficient. Obviously I did not express this clearly enough - sorry for that.
In the meantime I tried to determine how many loops it will take for a relative error below 5E-35. So I wrote a short routine for Casio's Keisan service and set it to 38 digits working precision (sorry, no Mathematica or Maple here). Finally I tried to describe the relation between x and the required number of terms n, and here's the result: n >= 12 + 940 / (x - 1/3)^1,57or alternatively n >= 12 + 1000 / (x - 0,3)^1,6In both cases this gives a number of terms slightly higher than required (in a perfect calculator ;-)). The formula works down to x = 1 so that it can even be used if the break point is set to 1,28155. Which is what I now use on the emulator - it is fast enough for 1000 terms in virtually no time. ;-) Edit: I also tried |x| > 3,14 with the emulator's function, and yes, Phi(-3,2) has only 27 correct digits (when rounded). The user code function, based on the CF method with the number of terms calculated with the formula above, agrees with the true value in all 34 digits. Well, OK, the last one is 1 while it actually should get rounded up to 2. ;-) I did some further tests and as far as I can tell, the error stays with a few units in the 34th digit. Until now, I did not notice any error exceeding 1 unit in the 33th digit. So if all this gets coded in C and 39 digit precision it looks like we can assume that all returned 34 digits should be exakt (within 1 ULP).
Dieter Edited: 10 May 2012, 1:54 p.m. after one or more responses were posted
Re: [WP34S] Curious Bug in Inverse Normal Function - Marcus von Cube, Germany - 05-10-2012 Dieter, I don't think we want to use exponentiation just to determine the number of loops so I implemented the old approach with 500/(|x|-1)+8 in DP. I reduced the numbers to 150 & 4 for single precision but I didn't do any tests.
Re: [WP34S] Curious Bug in Inverse Normal Function - Dieter - 05-10-2012 What? No exponents and powers ?-) n >= 300/(x-1,5) + 8This should be fine for 34 digits and x >= 2,326. It also does not overestimate the number of required loops as much as the 500/(x-1) method does. For 16-digit accuracy something like 100/(x-1) + 8 should be okay. But I have not tested this yet.
Edit: Here's my suggestion for 16-digit accuracy and x >= 2,326: n >= 85/(x-1,3) + 4 Finally: this does not have to be coded with high-precision DecNumber functions. Even standard C single precision will do. ;-) Dieter
Edited: 10 May 2012, 5:45 p.m. after one or more responses were posted
Re: [WP34S] Curious Bug in Inverse Normal Function - Marcus von Cube, Germany - 05-10-2012
Quote:They use LOGs which are are expensive. It's all about saving time, isn't it? I will try your formula but I cannot yet tell you when. Maybe Pauli is quicker. Can you try the latest incarnation and do some tests? It uses the 2.326 threshold. For SP, is 100/(x-1)+8 better than 150/(x-1)+4 ? Edited: 10 May 2012, 2:38 p.m.
Re: [WP34S] Curious Bug in Inverse Normal Function - Dieter - 05-10-2012 Quote:Sure. But an optimized n may save many useless loops. ;-)
Quote:The latter returns a loop count about 30-40% high. Replace the 150 by 110 and you get a very good estimate that works fine down to x = 1,7 and is just a few counts higher than the required minimum. I edited my last post to add another formula for 16-digit precision. According to my tests this should give no relative errors larger than 3E-17. Actually 110/(x-1) + 4 (and round up via CEIL) works just as well. #-)
You can try this yourself: Visit the Casio Keisan Website, go to the calculator and enter this: n = -int(-(110/(x-1) + 4));Press "Execute" and use "2 , 0.1 , 10" as your x-value. This will generate a table with the approximated and exact PDF as well as the relative error and finally the number of iterations. :-) I will also try the latest version in the emulator. For me this is always a lot of guesswork since I am not sure which files need to be updated. #-) My emulator version is 3.0 2725. Could you provide a list which files need to be replaced?
Dieter
Re: [WP34S] Curious Bug in Inverse Normal Function - Paul Dale - 05-10-2012 This is harder than you expect -- the routine doesn't know the requested accuracy. It is called from xrom so the machine will always be in double precision so that mode flag isn't helpful. The flag indicating the original user's single/double setting isn't adequate either since this code could be called from the QF routine that really does require some extra precision even though the user is running in single precision. We'll need an out of band way to indicate this. i.e. use one of the other xrom local flags to indicate what precision to use if the user's mode is set to single.
Re: [WP34S] Curious Bug in Inverse Normal Function - Marcus von Cube, Germany - 05-10-2012
Quote:Generally, wp34sgui.exe and wp34s-lib.dat from the windows bin directory should be fine. You may need to remove any older wp34s.dat file. If you want to use the assembler, update the tools directory contents, except for the .exe files which aren't up-to-date. You will need a working Perl installation. The best way is to install SVN, check out the complete tree once and then just do an SVN update to become up-to-date. EDIT: I've put up a new release file. This is the easiest option.
Edited: 11 May 2012, 4:45 a.m.
Re: [WP34S] Curious Bug in Inverse Normal Function - Marcus von Cube, Germany - 05-10-2012 You could just pass the required accuracy on the stack if necessary. We can always compute to an accuracy high enough even for nested calls in single precision. What is implemented now should be fine for all cases I can envision even if it may waste some cycles in single precision.
Edit: I've implemented Dieter's latest suggestions. Pauli, can you give an example when you need higher accuracy? This should be easy to implement. The nominator for SP is now 100 (instead of 83 as suggested by Dieter) which gives a little headroom. Edited: 11 May 2012, 4:48 a.m.
Re: [WP34S] Curious Bug in Inverse Normal Function - Paul Dale - 05-11-2012 Quote: I'm quite tempted to remove the special case for the standard normal distribution and just use gammap and gammaq for this. We'll lose speed in the quantile function but save lots of flash -- we lost almost three pages of flash special casing the normal cdf & pdf and doing the special purpose normal qf solver.
Re: [WP34S] Curious Bug in Inverse Normal Function - Paul Dale - 05-11-2012 More accuracy will be required when solving for the quantile function.
Re: [WP34S] Curious Bug in Inverse Normal Function - Dieter - 05-11-2012 I think we will need no extra headroom now. ;-) It's 3:45 a.m. in central Europe now - and I did some tests on the continued fraction approach for the Normal CDF and x > 2,326. This time I tried the Lentz algorithm in order to determine the exact number of terms required for a specific accuracy. This lead to a quite extensive table with exact data. And it turned out that the relation between x and the number of required terms can be described very precisely the following way, so that that we get the desired accuracy without significant overhead.
For single (16-digit) precision I think it's nice to have one or two guard digits to make sure that the last digit is rounded correctly. So I set the tolerance to 1E-18 to get 18 digits +/- 1 ULP. Which lead to... n >= 480 / x^2 + 9This matches almost exactly the number of required terms for the given accuracy. Only here and there the result is 1 or 2 terms lower, so that the relative error might increase to 1,2 E-18 or something similar.
For double (34-digit) precision I think it makes sense to simply use the full 39 internal digits. The last one or two may be slightly off anyway. With a target accuracy of 1E-39 I got this: n >= 2200 / x^2 + 20Just far out in the tails where only very few terms are required - and processed in virtually no time - the calculated n is a bit high due to the finally added constant. So I would definitely suggest that the number of terms for the Normal CDF should be set according to the respective formula above. They very exactly match the number of required terms for the desired accuracy while on the other hand no time is wasted for useless loops "just to be sure". ;-)
Dieter
Re: [WP34S] Curious Bug in Inverse Normal Function - Les Wright - 05-11-2012 Quote: For x > 0:
Lower Tail Normal CDF = 0.5*(1 + gammap(0.5*x^2, 0.5)) The parameter order is as you implement it in the calculator, and the opposite of the convention in Maple and Mathematica--the "a" parameter second, the "x" parameter first. See A&S 26.2.30.
Les
Re: [WP34S] Curious Bug in Inverse Normal Function - Paul Dale - 05-11-2012 Quote: I switched the parameter order for gammap and ibeta this morning :-)
Re: [WP34S] Curious Bug in Inverse Normal Function - Paul Dale - 05-11-2012 Quote:Lower Tail Normal CDF = 0.5*(1 + gammap(0.5*x^2, 0.5)) Thanks. I'll have a go with these two. I remember some difficulties with the lower tail cdf via the gammap function or rather via erf and then gammap. That's the way I did it originally.
Re: [WP34S] Curious Bug in Inverse Normal Function - Paul Dale - 05-11-2012 I remember the problem with the first formula. For x<<0 you can't just invert using 1 - cdf which cancellation. However, with both versions present, this isn't such an issue. - Pauli
Re: [WP34S] Curious Bug in Inverse Normal Function - Les Wright - 05-11-2012 Exactly. With x < 0, things are just reflected--the lower CDF formula for negative argument is just the upper CDF formula with positive statistic, and the upper tail formula is just the lower tail formula with positive statistic. You couldn't do that reflection thing when you didn't have gammaq available.
Les
Re: [WP34S] Curious Bug in Inverse Normal Function - Les Wright - 05-12-2012 I notice that in the documentation for the newly baptized P and Q the parameter order reflects what actually goes on the stack--second goes first so second is in y and the first is in x, just as the manual indicates. Works for me. With that in mind, in the normal distribution special cases x^2/2 goes into y and .5 into x before calling P or Q. So maybe I should've written things the conventional way to match the manual and the soon-to-be compiled functions.
Les
|