(50g) Normal Distribution

12032018, 04:20 PM
(This post was last modified: 01262019 10:04 PM by John Keith.)
Post: #1




(50g) Normal Distribution
I have been following these two threads and rather than trolling either with HP 50g programs I am starting this thread of normal distribution programs.
NOTE: The first two programs below should be considered obsolete. They are inaccurate for large values of x as pointed out by Dieter and Albert Chan in following posts. I am leaving them here for illustrative purposes only. The programs in posts 3 and 18 should be used instead. The function Q(x), the upper tail normal CDF, is implemented in the HP 50 as UTPN. It requires the mean on level 3, the variance on level 2, and x on level 1. Note that the HP Prime implements P(x), the lowertail CDF as NORMALD_CDF. Either function can be obtained from the other by subtracting the value from 1 as detailed in the above linked threads. For example, the lowertail normal CDF would simply (but inaccuratly) be Code:
The normal PDF Z(x), which is implemented on the Prime as NORMALD, is not built in to the HP 50. It can be easily (but inaccurately) calculated with the following program: Code:
UPDATE: Yes the normal PDF is included. It's called NDIST. Sorry to all for the confusion. Next, a program to generate normally distributed random numbers. This program uses the Marsaglia polar method and returns two random numbers: Code:
Finally, a somewhat useful example program that takes advantage of the fact that the above program returns two random numbers at a time. Given an integer or a list of array dimensions on level 1, the following program returns an array of normally distributed complex numbers. The program requires the above program stored in a variable named 'N2' in an accessible directory: Code:
I believe that all of the above programs will run on the 48g and 49g calculators but they have only been tested on the 50g. 

12032018, 07:05 PM
Post: #2




RE: (50g) Normal Distribution
(12032018 04:20 PM)John Keith Wrote: Either function can be obtained from the other by subtracting the value from 1 as detailed in the above linked threads. For example, the lowertail normal CDF would simply be Yes, that's simple, but not a good idea. Try this for x=–10 to see why. ;) A better method that returns accurate results is this one (in pseudocode – you will know how to do this on the 50g): Code: p = UPTN(abs(x)) By the way, another useful function is the twosided symmetric CDF A(x), i.e. the integral from –x to +x. Here similar considerations apply to get an accurate result for any x. If A(x) is calculated the tricky part is maintaining a decent accuracy for x close to 0, while for 1–A(x) the less trivial cases are those with large x. But this requires more than just a function for the upper or lower tail integral. What we need is an accurate value for the integral from 0 to x. Otherwise the trivial way of calculating A(x) = 1–2·UPTN(abs(x)) will suffer from digit cancellation for x<0,12566... As a rule of thumb, at x=10^{–n} the last n digits are lost. But in praxi this admittedly is not much of a problem. And for x<10^{–6} you can simply set A(x) = x·√(2/pi). ;–) (12032018 04:20 PM)John Keith Wrote: The normal PDF Z(x), which is implemented on the Prime as NORMALD, is not built in to the HP 50. It can be easily calculated with the following program: 3,14159265359? Does the 50g not have a constant for pi ?) Instead of 3,14159265359 you could directly code √(2pi) = 2,50662827463. But, more important, this direct method is not exact for large x. The reason for the problem is explained more detailled in the HP67/97 program thread you have linked to. If you take a look at the program you will see a more elaborate method for calculating Z(x) which provides better accuracy. (12032018 04:20 PM)John Keith Wrote: Next, a program to generate normally distributed random numbers. Does the 50g also offer the inverse Normal CDF, i.e. the quantile function? Then you can also use this function to generate normally distributed random numbers. Dieter 

12032018, 08:45 PM
(This post was last modified: 12032018 09:03 PM by John Keith.)
Post: #3




RE: (50g) Normal Distribution
Thank you for your valuable insights, Dieter.
(12032018 07:05 PM)Dieter Wrote: Yes, that's simple, but not a good idea. Try this for x=–10 to see why. ;) HP 50 version of above: Code:
(12032018 07:05 PM)Dieter Wrote: By the way, another useful function is the twosided symmetric CDF A(x), i.e. the integral from –x to +x. Here similar considerations apply to get an accurate result for any x. If A(x) is calculated the tricky part is maintaining a decent accuracy for x close to 0, while for 1–A(x) the less trivial cases are those with large x. But this requires more than just a function for the upper or lower tail integral. What we need is an accurate value for the integral from 0 to x. Otherwise the trivial way of calculating A(x) = 1–2·UPTN(abs(x)) will suffer from digit cancellation for x<0,12566... As a rule of thumb, at x=10^{–n} the last n digits are lost. But in praxi this admittedly is not much of a problem. And for x<10^{–6} you can simply set A(x) = x·√(2/pi). ;–) The HP 50 does have a constant pi but it is symbolic by default and converting it to a numeric value takes longer than the rest of the code to execute. \pi \>NUM saves about 5 bytes over the inline constant but I feel the extra memory is worth it if the code is to be executed inside a loop. Also, I cannot code √(2pi) directly because the actual formula is √(2*pi*v). The ... 2. * ... in the program is where the variance v is multiplied by 2. (12032018 07:05 PM)Dieter Wrote: But, more important, this direct method is not exact for large x. The reason for the problem is explained more detailled in the HP67/97 program thread you have linked to. If you take a look at the program you will see a more elaborate method for calculating Z(x) which provides better accuracy. The HP 50 does not have inverse CDF functions, nor does it have a builtin command for normal RANDs. The Prime does have all of the above functions but I am trying to provide the missing ones for the HP 50. I will look at your other suggestions when I have time but I fear some of them are "above my pay grade". John 

12032018, 09:53 PM
Post: #4




RE: (50g) Normal Distribution  
12032018, 10:18 PM
(This post was last modified: 12032018 10:21 PM by Dieter.)
Post: #5




RE: (50g) Normal Distribution
(12032018 09:53 PM)Albert Chan Wrote: Is it simpler just do this ? P(x) = Q(x) Yes, sure. Sometimes you don't realize the obvious. #) In this case there even is no noticeable difference in the accuracy of the result. I can't believe I didn't post this myself since I always use this in Excel spreadsheets with my own Q(x) code. #) Dieter 

12032018, 10:59 PM
Post: #6




RE: (50g) Normal Distribution  
12042018, 08:54 PM
(This post was last modified: 12042018 08:55 PM by John Keith.)
Post: #7




RE: (50g) Normal Distribution
So I have looked at the HP67 code for Z(x) (subroutine 'E') and I believe this is the RPL equivalent:
Code:
where the inline constant is SQRT(1/(2*pi)). I tried both the above code and my original prog with an input of 10 and the results differ by only one ULP. I am clearly not an expert on either statistics or the HP67 so I may be missing something here but I don't see an advantage for the more involved program above. 

12042018, 09:48 PM
(This post was last modified: 12042018 09:50 PM by Dieter.)
Post: #8




RE: (50g) Normal Distribution
(12042018 08:54 PM)John Keith Wrote: RPL equivalent: I don't know much about RPL, but 1/sqrt(2pi) is 0,39894228040143... However, on a 12digit calculator dividing by 2,50662827463 might be slightly more accurate. ;) (12042018 08:54 PM)John Keith Wrote: I tried both the above code and my original prog with an input of 10 and the results differ by only one ULP. I am clearly not an expert on either statistics or the HP67 so I may be missing something here but I don't see an advantage for the more involved program above. An input of x=10 will not make much of a difference as x has merely two significant digits so that x²/2 is easily evaluated exactly. But try a larger value with more digits, say 61/3 = 20,3333333333. The exact PDF for this (12digit rounded) input is 6,64644887122 E91. What do you get with the direct method? Dieter 

12042018, 11:06 PM
(This post was last modified: 12052018 02:04 AM by Albert Chan.)
Post: #9




RE: (50g) Normal Distribution
(12042018 09:48 PM)Dieter Wrote: ... try a larger value with more digits, say 61/3 = 20,3333333333. You really pick a hard one ... x²/2 = 206.722 222221 544444 444445 ~ 206.722 222222, almost a halfway case Doing everything with Casio FX115MS, which have 12 digits precision Direct method, with k = √(2 Pi): Z(x) = exp(x²/2) ÷ k = 6.64644 886816e91, error = 306 ULP 3 Exp method, I = INT(x) = 20, F = FRAC(x) = 0.33333 33333: Z(x) = exp(I²/2) exp(I F) exp(F²/2) ÷ k = 6.64644 887123e91, error = +1ULP 2 Exp method: B = INT(x²)/2 = 413/2 = 206.5 Z(x) = exp(B) exp(B  I²/2  I F  F²/2) ÷ k = 6.64644 887123e91, error = +1 ULP I think direct method is good enough (9 digits precision, even with extreme case) 

12052018, 03:02 PM
(This post was last modified: 12052018 04:42 PM by Albert Chan.)
Post: #10




RE: (50g) Normal Distribution
For x = 16.4285 714286 (~ 115/7), Z(x) errors is much worse.
x²/2 suffered from double roundings, 134.948 979592 30612 ... => 134.948 979593 2 Exp method got Z(x) correctly = 9.84720 298683 e60 Direct method shorted 687 ULP > 9.84720 297996 e60 Edit: above assumed calculator do halfroundup. With halfwaytoeven, direct method got 9.84720 298986 e60, error = +303 ULP 

12052018, 08:49 PM
(This post was last modified: 12052018 08:52 PM by Dieter.)
Post: #11




RE: (50g) Normal Distribution
(12042018 11:06 PM)Albert Chan Wrote: I think direct method is good enough (9 digits precision, even with extreme case) The point is not that you get nine digits, but that you lose three. The program in the HP67/97 thread with the improved was designed to achieve close to machine precision for most cases, i.e. roughly 9 out of 10 digits. If one intermediate result returns only 7 valid digits this is an absolute nogo. Dieter 

12062018, 02:58 PM
Post: #12




RE: (50g) Normal Distribution
(12052018 03:02 PM)Albert Chan Wrote: For x = 16.4285 714286 (~ 115/7), Z(x) errors is much worse. I have compared your 2 EXP method and Dieter's original 3 Exp method and the 3 Exp version is actually faster on the HP 50. Exp is pretty fast on the HP 50 and the 3 Exp program is shorter, which helps because of RPL system overhead. Both methods give results accurate to +/ 1 ULP for all cases stated above. The one problem I am having, due to my rusty algebra skills, is factoring in the variance. I would like the program to take the same arguments as the UTPN function, mean on level 3, variance on level 2, and x on level 1, as my attempt in the first post did. However, if I divide by 2v instead of by 2, I get the wrong answer. Any helpful hints would be appreciated. 

12062018, 08:30 PM
(This post was last modified: 12072018 11:52 PM by Albert Chan.)
Post: #13




RE: (50g) Normal Distribution
Hi, John Keith
You are right. On a calculator, 3 Exp method maybe faster. I did the test on the computer. Cost of PDF function is dominated by Exp. 2 Exp method take 2X as long, 3 Exp method take 2.83X (using C language) "Splitting" factor 1/√(2 pi) inside Exp reduce error somewhat: PHP Code: def pdf(x, k=0.08106146679532726, exp=math.exp, floor=math.floor): 

12062018, 10:54 PM
Post: #14




RE: (50g) Normal Distribution
(12062018 08:30 PM)Albert Chan Wrote: Using lrint() and move factor √(2 pi) inside Exp reduce error somewhat: What does this lrint() function do, and what programming language are we talking about in the first place? (12062018 08:30 PM)Albert Chan Wrote: const double k = 0.9189385332046728; /* ln(sqrt(2 pi)) */ That's a border case – ln(sqrt(2*pi)) does not round very well to 16 decimals. The exact value is 0,91893 85332 04672 74178... Maybe one has to find out which decimal value is closer when the binary representation is considered. ;) Dieter 

12062018, 11:08 PM
Post: #15




RE: (50g) Normal Distribution
(12062018 02:58 PM)John Keith Wrote: The one problem I am having, due to my rusty algebra skills, is factoring in the variance. I would like the program to take the same arguments as the UTPN function, mean on level 3, variance on level 2, and x on level 1, as my attempt in the first post did. However, if I divide by 2v instead of by 2, I get the wrong answer. Any helpful hints would be appreciated. What is "2v"? To get the Normal PDF for a distribution with mean µ and variance σ², normalize the variable x first by calculating (x–µ)/σ. Then apply the Standard Normal PDF Z(), and finally divide that value by σ. In other words: simply calculate 1/σ · Z[(x–µ)/σ]. Dieter 

12062018, 11:36 PM
(This post was last modified: 12072018 12:00 AM by Albert Chan.)
Post: #16




RE: (50g) Normal Distribution
Hi, Dieter
lrint() is from C language, roundtointeger. For rounding mode = FE_TONEAREST, it is same as nearestinteger, using bankers rounding. My latest update change the code to more familiar Python. (12062018 10:54 PM)Dieter Wrote:(12062018 08:30 PM)Albert Chan Wrote: const double k = 0.9189385332046728; /* ln(sqrt(2 pi)) */ I just switched to k = 1  ln(sqrt(2 pi)) instead. This constant is more precise, gaining 5 bits precision (or 1 decimal digit). Round to 16 digits, 0.08106 14667 95327 25821 ... ==> 0.08106 14667 95327 26 This happens to be same value for binary float 

12072018, 11:21 PM
(This post was last modified: 12082018 06:56 PM by John Keith.)
Post: #17




RE: (50g) Normal Distribution
(12062018 11:08 PM)Dieter Wrote: What is "2v"? By 2v I meant 2 * variance as in the formula I found in several references, and which was the basis of my program in the first post: 1/Sqrt(2*pi*v)*e^((xm)^2/(2*v)) (Please excuse my lack of MathJax ability.) I did manage to figure out the algebra well enough to get Albert's 2Exp method working properly with mean and variance. It is a big, slow beast of a program, though, and I will study your post carefully and hopefully come up with something better. Thanks for your help, John 

12082018, 07:13 PM
Post: #18




RE: (50g) Normal Distribution
Thanks to Albert and Dieter's help I have a working program for the normal PDF with mean and variance. Dieter's method for dealing with the mean and variance greatly simplified the program.
On the other hand, the idea of using 1ln(sqrt(2*pi)) proved to be less accurate on the HP 50 so I stuck with dividing by sqrt(2*pi) instead. The resulting program is still big and ugly but reasonably fast, around 43 milliseconds. It is accurate to +/ 1 ULP for all test cases in the above posts. Code:


12082018, 08:15 PM
(This post was last modified: 12082018 08:18 PM by Albert Chan.)
Post: #19




RE: (50g) Normal Distribution
Hi, John Keith
Now that you have a good PDF function, here is the bad news. The extra precision is basically useless. zscore = (x–µ)/σ, with all 3 variables nonexact. Even if you have accurate zscore, say 20.33333 33333 (12 digits precision) eps = ±5e11, thus PDF = 6.64644886446e91 to 6.64644887799e91 PDF value is good only for 9 digits, even if the function can generate more. That was why I mentioned direct method is good enough, 6.64644887e91 

12082018, 09:18 PM
Post: #20




RE: (50g) Normal Distribution
(12082018 08:15 PM)Albert Chan Wrote: The extra precision is basically useless. Are they? For the standard Normal distribution µ is exactly 0 and σ is exactly 1. So there is a real benefit if the calculation is done with the best possible accuracy. Sure, the general assumption for all numeric calculations, be it on a HP67, HP50g or a contemporary computer, is that the input is exact. So the input is not √2 but 1,41421356237, and not pi but 3,14159265359. That's a limitation that always applies. But in other cases the input can be exact. And here there improved PDF offers a benefit. Dieter 

« Next Oldest  Next Newest »

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