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: Yesterday 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: Yesterday 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. 

Yesterday, 10:49 PM
(This post was last modified: Yesterday 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. 

Today, 12:59 AM
(This post was last modified: Today 02:11 AM by Albert Chan.)
Post: #30




RE: Accurate Normal Distribution for the HP67/97
(Yesterday 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.43280 95260 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^{18} ~ 0.000147, y^{4}/24 < 2.0e17 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, match true value 

« Next Oldest  Next Newest »

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