Normal Probability Function

[CMarangon]

Hello!

I want to make a program to find P(x) and Q(x) for Normal Probability Function.
I found formulas for P(x) and Q(x) as shown in the picture below.
However they are hard integrals
and it seems that HP50 does not solve them.

I am looking for some kind of expanded formulas
like Taylor Polynomial.

Even Wolframalfa and Symbolab seem to be unable to help.

Can you help me?

[Joe Horn]

Not sure, but the 50g's built-in UTPN function seems to be a generalized version of P(x) in your table. The section of the 50g AUR about UTPN can be seen here: https://holyjoe.net/hp/UTPN.jpg

The 50g also has the NDIST(m,v,x) function which might be of some use to you. It returns the normal probability distribution (bell curve) at x based on the mean m and variance v.

Just guessing here... if not applicable, please ignore this reply.

[CMarangon]

(03-28-2022 09:11 PM)Joe Horn Wrote:  Not sure, but the 50g's built-in UTPN function seems to be a generalized version of P(x) in your table. The section of the 50g AUR about UTPN can be seen here: https://holyjoe.net/hp/UTPN.jpg

The 50g also has the NDIST(m,v,x) function which might be of some use to you. It returns the normal probability distribution (bell curve) at x based on the mean m and variance v.

Just guessing here... if not applicable, please ignore this reply.

Hi Joe!
I will take a look at UTPN and NDIST.
I was making a webpage for it.{
http://www.areaseg.com/gauss

I got these formulas from my Casio FX-5000F.
They work!

This is the corrected formula.
Thanks to Mike and Trojdor for help.
Ancient erros are in green.

[Mike T.]

This may be of some use https://www.sasview.org/docs/_downloads/.../sas_erf.c

For a translation into Visual BASIC take a look at the source code for the HP32E here https://www.hpmuseum.org/cgi-sys/cgiwrap...i?read=867 (specifically modStastics.bas)

[Paul Dale]

These functions are very difficult to compute accurately over their entire domain. The WP 34S uses the incomplete gamma function to work them out:

Lower Tail = 0.5*(1 + gammap(0.5*x^2, 0.5))
Upper Tail = 0.5*gammaq(0.5*x^2, 0.5)

This belies the difficulty computing gammap and gammaq accurately. The common approach is to break the domain into two pieces and compute the function using a series in one part and a continued fraction for the other. See e.g. section 6.5 in Abramowitz and Stegun. This fails for large arguments from memory and more recent approaches break the domain into more pieces. All this relies on being able to compute the (standard) gamma function accurately which is another problematic one -- albeit nowhere near as nasty as the imcomplete gamma function.

If you do implement the incomplte gamma function, you also get Poisson and χ² distributions for free.


Pauli

[Albert Chan]

(03-29-2022 02:53 AM)CMarangon Wrote:  I got these formulas from my Casio FX-5000F.
They work!

Yes ... but not well.
It seems formula converted from https://www.johndcook.com/blog/cpp_erf/
For upper tail area, to be more accurate, we use erfc(x) = 1 - erf(x)

ERFC(x) := horner([1.061405429,-1.453152027,1.421413741,-0.284496736,0.254829592,0],
1/(1+0.3275911*x)) * exp(-x*x)

P(x) := ERFC(x/sqrt(2))/2

However, formula is sensitive with constants, slight rounding may mess up everything.

With above setup, max absolute error about 7e-8
With Casio setup, max absolute error goes up to 83e-8

Worse case, at x = 0.32:
ERFC  P(x) = 0.37448 42341
Casio P(x) = 0.37448 49931
True  P(x) = 0.37448 41653

Error Ratio = (41653-49931) / (41653-42341) = -8278/-688 ≈ 12.

[Albert Chan]

(03-29-2022 12:00 PM)Paul Dale Wrote:  Upper Tail = 0.5*gammaq(0.5*x^2, 0.5)

Correction, factor off by 1/√(pi)

Γ(1/2, x²) = √(pi) * erfc(x)      // Abramowitz and Stegun, 6.5.17

P(x) = erfc(x/√2)/2 = Γ(1/2, x²/2) / (2*√(pi))

>>> from mpmath import *
>>> P1 = lambda x: erfc(x/sqrt(2))/2
>>> P2 = lambda x: gammainc(1/2,x*x/2) / (2*sqrt(pi))
>>> z = 1
>>> print P1(z), P2(z) # confirm P1 = P2
0.158655253931457 0.158655253931457

Quote:The common approach is to break the domain into two pieces and compute the function
using a series in one part and a continued fraction for the other.

For Xcas implementation of erf, see https://www.hpmuseum.org/forum/thread-15...#pid135589

[KeithB]

I love it when someone references "Big Red". 8^)

[CMarangon]

Hello everybody!

I made this program for HP50.
It works and is shorter than the one I made using the formula from Casio FX-5000F




Soon I will send it to HPCALC.

[CMarangon]

(03-29-2022 09:51 AM)Mike T. Wrote:  This may be of some use https://www.sasview.org/docs/_downloads/.../sas_erf.c

For a translation into Visual BASIC take a look at the source code for the HP32E here https://www.hpmuseum.org/cgi-sys/cgiwrap...i?read=867 (specifically modStastics.bas)

Hello!

I will use it in my PHP page to calculate G(x), P(x) and Q(x):-)

It has more numbers that the one I have got from my Casio calc.:-)

This is the page I made, unfortunbately HP50 functions like NDIST and others do not work in PHP

http://www.areaseg.com/gauss

I will also share the webpage code,
in few days.

Soon I ill translate it to English, but, for now, even in Potuguese, I believe that you all can understand.

[CMarangon]

hello!

Here is the source code for a PHP page

http://www.areaseg.com/gauss/gauss.zip

[Paul Dale]

As soon as any of these kinds of formulæ involve a 1 - term, you need to reconsider your choices. Accuracy will be lost when term approaches unity.

Pauli

[trojdor]

The primary formula listed has it's roots in Abramowitz and Stegun's
Handbook of Mathematical Functions, Tenth Printing p.932
Probability Functions 26.2.17

A modified version of this CDF approximation is used in many function libraries with good results.
(I use a similar modified approximation in my CDF programs for the HP-35s.)

However, your listing shows several typos in the constants, and that will adversely affect accuracy:

Denominator of f(t);
0.2316419 <--- your listing missing the final 9

For P(x), the constants are;
b1= 0.319381530
b2=-0.356563782
b3= 1.781477937 <--- your listing missing the extra 7 in the middle
b4=-1.821255978
b5= 1.330274429 <--- your listing has typo ... 3 as next to last digit should be 2

Additionally, modern CDF function libraries tend to replace the Sqrt(2*PI) in P(x) denominator with a pre-calculated constant, 2.50662827463, for speed / efficiency.

mike

[CMarangon]

Hi Mike and trojdor!

I will fix these mumbers.
Thank you!

[Albert Chan]

(03-30-2022 08:02 AM)trojdor Wrote:  Additionally, modern CDF function libraries tend to replace the Sqrt(2*PI) in P(x) denominator with a pre-calculated constant, 2.50662827463, for speed / efficiency.

It may be better to remove factor altogether, using P(x) = erfc(x/√(2))/2

Below formula practically matched erfc (A&S 7.1.26) based P(x) formula.
In other words, same max abs error of 7e-8.

P(x) := exp(-x*x/2) *
horner([0.53070271,-0.72657601,0.71070687,-0.14224837,0.1274148,0.], 1/(1+0.23164189*x))

Note: equation constants only required 8 digits precision
I discovered this using numpy to move "/2" into the coefs.

>>> a = [1.061405429,-1.453152027,1.421413741,-0.284496736,0.254829592,0]
>>> numpy.array(a) / 2
array([ 0.53070271, -0.72657601, 0.71070687, -0.14224837, 0.1274148 , 0. ])

Turns out, 8 digits is enough to get practically the same error plots.

[CMarangon]

Corrected Equation

You told it was wrong and I am not stubborn :-)

Now it has only accuracy errors, because of formula itself.

Thanks to Joe, Mike and "trojdor"!

Corrections are in green:

[Danito-Danone]

Hi all, From de HP-25 Applications Programs, here:
Handbook of Mathematical Functions With Formulas, Graphs,...
Edited by: Milton Abramowitz and Irene A. Stegun
National Bureau of Standards, 1968 / Dec. 1972 (Available on Internet)
Normal Distribution: Item 26.2.17 - page: 932
Inverse Normal integral: Item 26.2.23 - page 933
D-D