Accurate Normal Distribution for the HP67/97

12032018, 07:02 PM
Post: #21




RE: Accurate Normal Distribution for the HP67/97
(12032018 12:44 AM)Albert Chan Wrote: That is a cool formula ! I was expecting something much worse. Something much worse would not have fit into the limited HP67/97 memory. ;) (12032018 12:44 AM)Albert Chan Wrote: To prove above correction work: Ah, thank you very much. I have never tried a proof for that formula myself. But I wonder what that neglected fourth term with t^{4} will look like. Due to the limited hardware precision it will not make much sense, but anyway. I have a faint idea, but you will know better – Albert? BTW, the tricky part is calculating t in the first place. Especially evaluating Q(x)–p without too much digit cancellation. Dieter 

12032018, 08:25 PM
Post: #22




RE: Accurate Normal Distribution for the HP67/97
(12032018 12:44 AM)Albert Chan Wrote: Sum it all, and ignore O(h^{4}) terms, we get: I am too lazy to do 4th order correction by hand. Using Mathematica: z = c Exp[x²/2]; D[z, {x, 3}] / z ==> x (3 + x²) t = h  (x/2) h^2 + (x²1)/6 h^3  x*(x²3)/24 h^4; mess = Expand[t + (x/2) t^2 + (2x²+1)/6 t^3 + k t^4]; Simplify[(mess /. h^4 > c /. h > 0) / c] ==> k  1/24 x(7 + 6x²) 4th order x correction = t + x/2 · t² + (2x²+1)/6 · t³ + x*(6x²+7)/24 t^{4} 

12032018, 08:45 PM
(This post was last modified: 12032018 08:50 PM by Dieter.)
Post: #23




RE: Accurate Normal Distribution for the HP67/97
(12032018 08:25 PM)Albert Chan Wrote: 4th order x correction = t + x/2 · t² + (2x²+1)/6 · t³ + x*(6x²+7)/24 t^{4} Thank you very much. I just found an older VBA function in the Excel sheet where the initial approxmation for the Normal quantile was tested: Code: x = x + t * t * t * t * (6 * x ^ 3 + 7 * x) / 24 + t * t * t * (2 * x * x + 1) / 6 + t * t * x / 2 + t But I can't say how or when I got to that conclusion. ;) With four terms even a not too exciting fourdigit approximation like Hastings' yields 15 digit double precision accuracy with just one correction step. Probably even more  but Excel is limited to these 15 digits. Dieter 

12042018, 01:38 PM
(This post was last modified: 12042018 09:02 PM by Albert Chan.)
Post: #24




RE: Accurate Normal Distribution for the HP67/97
(06262016 09:34 PM)Dieter Wrote: Evaluating the PDF seems trivial, but accuracy may degrade significantly for large arguments of the exponential function. For example e^{1000/7} = 9,076766360 E63, but the 67/97 returns 9,076765971 E63. The error is caused by the fact that the fractional part of the argument carries only seven decimals which leads to an accuracy of merely seven significant digits. That's why the PDF is evaluated in a different way that requires three calls of e^{x}, but achieves better accuracy. It might be possible to reduce 3 calls of exp(x) to 2 Let B = INT(x²)/2, I = INT(x), F = FRAC(x): x²/2 = B + (B + I²/2 + I F + F²/2) Example: 10.23456789² / 2 = 52 + (52 + 10²/2 + 10*0.23456789 + 0.23456789²/2) ~ 52 + (52 + 50 + 2.3456789 + 0.02751104751) ~ 52 + 0.3731899475 exp(10.23456789² / 2) ~ exp(52) exp(0.3731899475) = 1.797267025e23 

12092018, 03:28 PM
(This post was last modified: 12122018 04:43 PM by Albert Chan.)
Post: #25




RE: Accurate Normal Distribution for the HP67/97
We can use Taylor expansion to calculate Z(x+h) with 1 Exp call
Z(x+h) = Z(x) + Z'(x) h + Z''(x)/2 h² + Z'''(x)/6 h³ + ... Z(x+h) / Z(x) = 1  x h + (x²  1)/2 h²  x (x²  3)/6 h³ + ... Example: calculate Z(20.3333 333333) x = zscore rounded to FIX4, thus exp(x²/2) can be calculate accurately. h = zscore  x, thus h ≤ 5e5 > above example, x=20.3333, h=3.33333e5 > Z(x) = exp(x²/2) / √(2 Pi) = 6.65095 520534 e91 Correction to get to Z(x+h): Code: c1 = x h = 677775.98889 e9 Z(x+h) ~ Z(x) + c Z(x) = 6.64644 887122 e91, matched true value 

12092018, 08:17 PM
(This post was last modified: 12102018 10:25 AM by Albert Chan.)
Post: #26




RE: Accurate Normal Distribution for the HP67/97
Above 1 Exp method, if done on the computer, is almost as fast as direct method.
In other words, the extra precision almost without cost. PHP Code: double pdf(double z) // 1 Exp method Updated: fix typo in constant 1/√(2 Pi) 

12092018, 08:47 PM
Post: #27




RE: Accurate Normal Distribution for the HP67/97  
12112018, 02:33 PM
(This post was last modified: 12122018 02:52 PM by Albert Chan.)
Post: #28




RE: Accurate Normal Distribution for the HP67/97
Discovered a simpler and faster 1 Exp method.
Instead of correcting Z(x) to Z(x+h) via Taylor expansion of Z(x+h), do this: Z(x+h) / Z(x) = exp(½(x+h)² + ½(x²)) = exp(x h  h²/2) = exp(y) exp(y) = 1 + y + y²/2 + y³/6 + ... y ~ x h, a very small number, above ... can be ignored. PHP Code: double pdf(double z) // 1 Exp method, revised Note: for calculator program, return z * (1. + y) instead Calculator might underflow the correction term. 

12122018, 10:49 PM
(This post was last modified: 12122018 10:50 PM by John Keith.)
Post: #29




RE: Accurate Normal Distribution for the HP67/97
That is a neat program, although I can only begin to understand it. I tried it on the HP 50 and while it is a smaller and faster program than the 2Exp program in this thread, it is not as accurate if the mean and variance are other than 0 and 1 respectively. For "standard" values of x it is accurate to +/ 1 ULP.
My implementation is straightforward except for the second line, double x = (float) z; which turned into 6 RND i.e. round to 6 digits. I did not experience any underflow, just a less accurate result. 

12132018, 12:59 AM
(This post was last modified: 12142018 09:55 PM by Albert Chan.)
Post: #30




RE: Accurate Normal Distribution for the HP67/97
(12122018 10:49 PM)John Keith Wrote: My implementation is straightforward except for the second line, Round to 6 digits does not work. Example: 4.56789²/2 = 10.4328 095260 5 Try round z to FIX4, we get: y < 21.4 * 5e5 = 0.00107, y^{4}/24 < 5.5e14 3rd order correction also enough for C double: From header file float.h: DBL_EPSILON = 2^{52} ~ 2.22e16 y < 38.6 * 2^{19} ~ 0.000074, y^{4}/24 < 1.2e18 Quote:I did not experience any underflow, just a less accurate result. An example of underflow, try Z(20.5000 00002): z = Z(20.5) = 2.21198 438021 e92 y = 0.5*(20.5  20.5000 00002)*(20.5 + 20.5000 00002) ~ 4.1e8 Other correction terms y²/2, y³/6 are too small, and can be ignored. > correction = y z ~ 9.069e100, underflow to 0 That was why I mentioned do z*(1. + y) for calculator program. > Z(20.5000 00002) = Z(20.5) * 0.99999 9959 = 2.21198 428952 e92, matched true value 

12132018, 08:43 PM
(This post was last modified: 12132018 09:06 PM by Dieter.)
Post: #31




RE: Accurate Normal Distribution for the HP67/97
(12132018 12:59 AM)Albert Chan Wrote:(12122018 10:49 PM)John Keith Wrote: My implementation is straightforward except for the second line, That's not what John is doing. Look closer: it's not 6 RND but –6 RND. Negative arguments round to a certain number of significant digits! It's like SCI 5 RND on a classic RPN device. Edit: I just noticed that your example also rounds to 6 significant digits. So you are correct here. If the division by 2 also has to be exact, the argument should be rounded to 5 significant digits instead, i.e. –5 RND. I think this is better than FIX 4 as the latter again rounds to 6 significant digits for all x≥10. BTW the largest x is not 38,6 but 47,918. The working range of the 12digit HPs is larger than IEEE 754 double precision, it goes down to 1E–499. You will know what this means for the question we're discussing here. ;) Dieter 

12132018, 10:02 PM
Post: #32




RE: Accurate Normal Distribution for the HP67/97
(12032018 03:53 AM)Valentin Albillo Wrote:(12022018 07:44 PM)Albert Chan Wrote: I do not have access to the Abramovitz & Stegun book. As valentin says. Also wiki https://en.wikipedia.org/wiki/Abramowitz_and_Stegun "Because the Handbook is the work of U.S. federal government employees acting in their official capacity, it is not protected by copyright in the United States. While it could be ordered from the Government Printing Office, it has also been reprinted by commercial publishers, most notably Dover Publications (ISBN 0486612724), and can be legally viewed on and downloaded from the web. " Some torrent with files related to calculators have it as well. Wikis are great, Contribute :) 

12132018, 10:53 PM
Post: #33




RE: Accurate Normal Distribution for the HP67/97
(12132018 10:02 PM)pier4r Wrote:(12032018 03:53 AM)Valentin Albillo Wrote: You can download it for free using this link (470page PDF, perfectly legal download): For the record: the 470 page version is a really nice and clean scan, but compared to the printed book it is a heavily stripped down version that lacks most of the tables. Yes, nowadays most of these are obsolete, but they add to the somewhat nostalgic feelings you may have while you look at the pages of this book. ;) That's why I prefer a different PDF version. With 53 MB It's not much larger than the 470page version, but it offers no less than 1060 pages (!). No, I can't remember from where I got this one. But I really like all these tables, e.g. table 23.4 with the sum of positive powers Σk^{n}. This is on page 817 which is missing in the stripped down version. Did you know that the sum of the 10th powers of 1...100 equals 9 59924 14243 42419 24250 ?) Dieter 

12132018, 11:13 PM
Post: #34




RE: Accurate Normal Distribution for the HP67/97
Hi, Dieter
I did not know about the ±499 exponent range. That is huge ! FIX4 was based on ±99 exponent range, we have Z(x) = 0.0 if x ≥ 21.4 So, FIX4 work: 5 RND for 1 ≤ x < 10, 6 RND for x ≥ 10 With ±499 exponent range, FIX4 is too much. Luckily, with revised 1 Exp method, adding correction terms is easy.


12142018, 02:43 AM
(This post was last modified: 12192018 01:39 AM by Albert Chan.)
Post: #35




RE: Accurate Normal Distribution for the HP67/97
Another 1 Exp method, but correct from direct method: (Revision 2)
Example: Z(z = 20.3333 333333), using Emu48 B = z²/2 = 206.722 222222 D = exp(B) / √(2 Pi) = 6.64644 886819 e91, error = 303 ULP D is very close to true result, only linear correction is enough ! x = z rounded to 5 digits = 20.333 h = z  x = 0.0003 333333 Z(z) = D (B  x²/2  x h  h²/2 + 1) = D * 1.00000 000046 = 6.64644 887125 e91, error = +3 ULP Edit: if correction underflow is not an issue, D + D (B  x²/2  x h  h²/2) may be more accurate. We get D + D * 45.5555 555 e11 = 6.64644 887122 e91, matched true result. 

12142018, 07:42 PM
(This post was last modified: 12142018 07:42 PM by Dieter.)
Post: #36




RE: Accurate Normal Distribution for the HP67/97  
12142018, 08:36 PM
Post: #37




RE: Accurate Normal Distribution for the HP67/97
(12142018 07:42 PM)Dieter Wrote: The WP34s (in DP mode) and Free42 Decimal support numbers as low as 1E–6143. ;) Though it's hard to imagine any application where the number of significant digits matter when using numbers of this magnitude... [preparing to be educated how this may be naive...] Bob Prosperi 

12142018, 10:25 PM
Post: #38




RE: Accurate Normal Distribution for the HP67/97
(12132018 08:43 PM)Dieter Wrote: That's not what John is doing. Look closer: it's not 6 RND but –6 RND. I did try both 5 RND and 7 RND and both made the error much worse. I did not try other values, however. I just saw Albert's latest posts and haven't tried any of the ideas there. 

12152018, 12:06 PM
(This post was last modified: 12152018 12:07 PM by pier4r.)
Post: #39




RE: Accurate Normal Distribution for the HP67/97
(12132018 10:53 PM)Dieter Wrote: Yes, nowadays most of these are obsolete, but they add to the somewhat nostalgic feelings you may have while you look at the pages of this book. ;) I don't think they are useless (and I remember my primary school math books had them, then they disappeared). Why do I think that math tables are somewhat useful today? It is a bit like paper dictionaries. You search a word and you may stumble in others while you search it. Or you can do a bit of random looking. Of course digital dictionaries can simulate this, but not as good as paper ones yet. The same with the mathematical tables. Sometimes curiosity arises from what one experiences. A clear table with certain numbers may induce the reader to explore more. With a calculator, yes, one can reproduce all the tables one wants (actually it could be quite a relaxing task) but one doesn't think about a result until one needs it. Therefore one doesn't have the table of numbers already compiled in front of him to start to notice patterns or ask questions. Aside from the above, I believe one could find mathematical tables in PDF or, as I said, one can make a little relaxing project to compile them with a calculator (not a computer!). Especially for large numbers it may be challenging. Actually it won't be much relaxing if one wants to be sure there are no typos or errors. I remember Jürgen posted a video where a 19C printed the first 10'000 prime numbers (n1). Well one can go and print (maybe in a digital way?) the results of the mathematical tables. n1: http://www.hpmuseum.org/forum/thread10243.html Wikis are great, Contribute :) 

12152018, 06:38 PM
(This post was last modified: 12152018 06:44 PM by Dieter.)
Post: #40




RE: Accurate Normal Distribution for the HP67/97
(12142018 10:25 PM)John Keith Wrote: I did try both 5 RND and 7 RND and both made the error much worse. I did not try other values, however. I just saw Albert's latest posts and haven't tried any of the ideas there. Then you have found your best way to round the input. And indeed rounding to 6 significant digits may be the best choice here. The square of such a number is exact, and the division by 2 will still be exact if the mantissa of x² does not exceed 2. This is true if the mantissa of x does not exceed √20, i.e. up to 4,47213. Which accounts for 44,72% of all possible x. The remaining 55,28% are also exact if the final digit of x² is even. Otherwise the value of x²/2 may be 1/2 ULP high. So an error only occurs in only 27,64% of all cases, and even if this happens the result if only half an ULP off. I have not analyzed what this means for the final Z(z) result, but a slight error in the last digit may always remain in such calculations. After all there is no guarantee that the e^{x} function is accurate to half an ULP. Take a look at the 15C Advanced Functions Handbook where the possible errors of different functions are explained in detail. Here an error of slightly more than 1 ULP (but less than 3 ULP) is possible. Please correct me if I'm wrong here. ;) Dieter 

« Next Oldest  Next Newest »

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