(50g) Normal Distribution

[John Keith]

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 lower-tail 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 lower-tail normal CDF would simply (but inaccuratly) be

Code:


\<< UTPN 1. SWAP - \>>

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:


\<< UNROT 2. * UNROT - SQ OVER / NEG EXP SWAP 3.14159265359 * \v/ /
\>>

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:


\<<
  WHILE RAND 2. * 1. - RAND 2. * 1. - DUP2 SQ SWAP SQ + DUP 1. \>=
  REPEAT DROP2 DROP
  END DUP LN -2. * SWAP / \v/ SWAP OVER * UNROT *
\>>

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:


\<< DUPDUP \-> d
  \<< TYPE 5. SAME { 1. + \PILIST } IFT # 1d SWAP R\->B
    START N2 R\->C
    NEXT d \->ARRY
  \>>
\>>

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.

[Dieter]

(12-03-2018 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 lower-tail normal CDF would simply be

Code:

\<< UTPN 1. SWAP - \>>

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))
if x>0 then p = 1-p
return p

By the way, another useful function is the two-sided 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). ;–)

(12-03-2018 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:

Code:

\<< UNROT 2. * UNROT - SQ OVER / NEG EXP SWAP 3.14159265359 * \v/ /
\>>

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.

(12-03-2018 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

[John Keith]

Thank you for your valuable insights, Dieter.

(12-03-2018 07:05 PM)Dieter Wrote:  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))
if x>0 then p = 1-p
return p

HP 50 version of above:

Code:


\<< UNROT PICK3 ABS UTPN SWAP 0. > { 1. SWAP - } IFT
\>>

(12-03-2018 07:05 PM)Dieter Wrote:  By the way, another useful function is the two-sided 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). ;–)



3,14159265359? Does the 50g not have a constant for pi ?-)
Instead of 3,14159265359 you could directly code √(2pi) = 2,50662827463.

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 in-line 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.

(12-03-2018 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.



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

The HP 50 does not have inverse CDF functions, nor does it have a built-in 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

[Albert Chan]

(12-03-2018 07:05 PM)Dieter Wrote:  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))
if x>0 then p = 1-p
return p

Is it simpler just do this ? P(x) = Q(-x)

[Dieter]

(12-03-2018 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

[John Keith]

(12-03-2018 09:53 PM)Albert Chan Wrote:  Is it simpler just do this ? P(x) = Q(-x)

Yes, in theory. However, Q(10) is about 7.62E-24, and the HP 50 returns 1. for Q(-10). This is to be expected since 1 - 7.62E-24 is 1 to machine precision, but not very informative.

[John Keith]

So I have looked at the HP-67 code for Z(x) (subroutine 'E') and I believe this is the RPL equivalent:

Code:


\<< DUP IP SQ 2. / NEG EXP OVER IP PICK3 FP * EXP / SWAP FP SQ EXP \v/ / .3989422804028 *
\>>

where the in-line 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 HP-67 so I may be missing something here but I don't see an advantage for the more involved program above.

[Dieter]

(12-04-2018 08:54 PM)John Keith Wrote:  RPL equivalent:

Code:

\<< DUP IP SQ 2. / NEG EXP OVER IP PICK3 FP * EXP / SWAP FP SQ EXP \v/ / .3989422804028 *
\>>

where the in-line constant is SQRT(1/(2*pi)).

I don't know much about RPL, but 1/sqrt(2pi) is 0,39894228040143...
However, on a 12-digit calculator dividing by 2,50662827463 might be slightly more accurate. ;-)

(12-04-2018 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 HP-67 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 (12-digit rounded) input is 6,64644887122 E-91. What do you get with the direct method?

Dieter

[Albert Chan]

(12-04-2018 09:48 PM)Dieter Wrote:  ... try a larger value with more digits, say 61/3 = 20,3333333333.
The exact PDF for this (12-digit rounded) input is 6,64644887122 E-91.

You really pick a hard one ...
x²/2 = 206.722 222221 544444 444445 ~ 206.722 222222, almost a half-way case

Doing everything with Casio FX-115MS, which have 12 digits precision

Direct method, with k = √(2 Pi):

Z(x) = exp(-x²/2) ÷ k = 6.64644 886816e-91, 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 887123e-91, 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 887123e-91, error = +1 ULP

I think direct method is good enough (9 digits precision, even with extreme case)

[Albert Chan]

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 e-60
Direct method shorted 687 ULP -> 9.84720 297996 e-60

Edit: above assumed calculator do half-round-up.
With halfway-to-even, direct method got 9.84720 298986 e-60, error = +303 ULP

[Dieter]

(12-04-2018 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 no-go.

Dieter

[John Keith]

(12-05-2018 03:02 PM)Albert Chan Wrote:  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 e-60
Direct method shorted 687 ULP -> 9.84720 297996 e-60

Edit: above assumed calculator do half-round-up.
With halfway-to-even, direct method got 9.84720 298986 e-60, error = +303 ULP

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.

[Albert Chan]

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):
    b = 0.5 * floor(x * x) # k = 1 - ln(√(2 pi))
    xi = floor(x)          # same k for binary or decimal float
    xf = x - xi
    return exp(-b - 1.) * exp(b - 0.5*xi*xi - xi*xf - 0.5*xf*xf + k) 


[Dieter]

(12-06-2018 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?

(12-06-2018 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

[Dieter]

(12-06-2018 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

[Albert Chan]

Hi, Dieter

lrint() is from C language, round-to-integer.
For rounding mode = FE_TONEAREST, it is same as nearest-integer, using bankers rounding.
My latest update change the code to more familiar Python.

(12-06-2018 10:54 PM)Dieter Wrote:  

(12-06-2018 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. ;-)

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

[John Keith]

(12-06-2018 11:08 PM)Dieter Wrote:  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

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^-((x-m)^2/(2*v))

(Please excuse my lack of MathJax ability.)

I did manage to figure out the algebra well enough to get Albert's 2-Exp 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

[John Keith]

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 1-ln(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:


\<< ROT - SWAP \v/ SWAP OVER / DUP SQ IP .5 * UNROT DUP IP SWAP FP \-> s i f
  \<< DUP NEG EXP SWAP i SQ .5 * - i f * - f SQ .5 * - EXP * 2.50662827463 s * /
  \>>
\>>

[Albert Chan]

Hi, John Keith

Now that you have a good PDF function, here is the bad news.
The extra precision is basically useless.

z-score = (x–µ)/σ, with all 3 variables non-exact.

Even if you have accurate z-score, say 20.33333 33333 (12 digits precision)

eps = ±5e-11, thus PDF = 6.64644886446e-91 to 6.64644887799e-91

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.64644887e-91

[Dieter]

(12-08-2018 08:15 PM)Albert Chan Wrote:  The extra precision is basically useless.

z-score = (x–µ)/σ, with all 3 variables non-exact.

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