+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Software Libraries (/forum-10.html)
+--- Forum: General Software Library (/forum-13.html)
+--- Thread: (50g) Normal Distribution (/thread-11898.html)
Pages: 1 2
RE: (50g) Normal Distribution - John Keith - 12-08-2018 09:28 PM
(12-08-2018 08:15 PM)Albert Chan Wrote: Hi, John Keith
Now that you have a good PDF function, here is the bad news.
The extra precision is basically useless.
Thanks, Albert. I realize that 12 digits of precision are overkill for any real-world application. After all, if you know that the chance of some outcome is 1 in 1E60 it hardly seems necessary to know it to 12 digits! Mainly, I just like the idea of a user function being as accurate as the built-in commands.
I have also learned a lot about efficient programming and numerical accuracy from your and Dieter's posts, and those of many others as well.
RE: (50g) Normal Distribution - Albert Chan - 12-15-2018 03:03 PM
(12-08-2018 07:13 PM)John Keith Wrote: ... the idea of using 1-ln(sqrt(2*pi)) proved to be less accurate on the HP 50 so I stuck with dividing by sqrt(2*pi) instead.
For 12-digits calculator, this constant is better, relative error ~ 1.25e-18
k = e-0.88453 / √(2 Pi) = 0.164726 536723 000000 206 ... => 0.164726 536723
Z(x) = exp(0.88453 - B) exp(B - I²/2 - I F - F²/2) * k
Examples:
Z(20.33333 33333) = exp(0.88453 - 413/2) exp(-0.222222 221544) * 0.164726 536723 = 6.64644 887123 e-91
Z(16.4285 714286) = exp(0.88453 - 269/2) exp(-0.448979 592306) * 0.164726 536723 = 9.84720 298683 e-60
edit: 1 - ln(sqrt(2*pi)) constant had a flaw I missed. exp(-1-B) term might underflow to 0.0
RE: (50g) Normal Distribution - Albert Chan - 12-15-2018 08:26 PM
(12-15-2018 03:03 PM)Albert Chan Wrote: k = e-0.88453 / √(2 Pi) = 0.164726 536723 000000 206 ... => 0.164726 536723
Combine above k with 1 Exp method, revision 2, calculate Z(z):
B = z*z/2
D = exp(0.88453 - B) * 0.164726 536723
x = z, rounded to 5 digits
h = z - x
y = B - x²/2 - x h - h²/2
Z(z) = D + y D
Examples:
Z(20.33333 33333) = 6.64644 886819 e-91 + 303 ULP = 6.64644 887122 e-91
Z(16.4285 714286) = 9.84720 298984 e-60 − 301 ULP = 9.84720 298683 e-60
RE: (50g) Normal Distribution - John Keith - 01-26-2019 10:01 PM
Update: It turns out I was wrong, the HP 50g (48G and 49 as well) does have a command for the normal PDF. It's called NDIST and its right there in the MATH/PROBABILITY menu, 2nd page, along with the other distribution-related commands. DOH!
Still, we all learned a lot about numerical accuracy in this thread, right?
RE: (50g) Normal Distribution - pier4r - 01-26-2019 10:12 PM
I think sharing alternative ways, even if solutions already exists, is always helpful.
- for discussions (see precision and so on)
- to collect info that solutions already exist (see ndist) that may not be known to others.
- to show examples how a problem can be solved, approach that is not easy to see in a compiled library (or even a library with difficult to read open source code)
In the minimal case, at least the person sharing the "redundant" solution is now part of the club "I did it too". Same as the club of "I solved it too" when there is a problem to solve of which the solution is known.
Therefore kudos.