(15C)(DM15) - PDF/CDF and Inverse of a Normal Distribution

[deetee]

Hello!

This is my second program I wrote for my new DM15.

It calculates the probability density function (PDF), the cumulative distribution function (CDF) and their inverse functions (!) of a standardized normal function.

Call it with the argument (Y) and the desired function number (X). The function number is 1 for the PDF, 2 for its inverse, 3 for the CDF and 4 for its inverse.

To calculate the CDF I used some approximating function rather than integrating the PDF, which offers a accuracy of at least 2 digits (error < 1%) and is much faster.

The label to call is D (distribution).

Examples:
* 0 (Y) and 1 (X) calculates PDF(0)=0.39
* 0.1 (Y) and 2 (X) calculates 1.66 (the positive solution) where the PDF is 0.1
* 1 (Y) and 3 (X) calculates CDF(1)=0.84
* 0.95 (Y) and 4 (X) calculates 1.64 where CDF(1.64)=0.95

The Code:

Code:


# ------------------------------------------------------------------------------
#T:Normal Distribution
#D:Calculates PDF, PDFinverse, CDF and CDFinverse of a standardized normal distribution.
#D:
#D:INPUT:
#D:Y: Argument
#D:X: Function
#D:
#D:FUNCTION:
#D:1: PDF
#D:2: PDFinverse, Argument = 0 ... 1/sqrt(2*pi)
#D:3: CDF
#D:4: CDFinverse, Argument = 0 ...1
#D:
#D:PDF(x) = 1/sqrt(2*PI) * exp(-x*x/2)
#D:
#D:CDF(x) = (integral, [-inf;x]) PDF(z) * dz =(approx)= 1/(1 + exp(-0.07*x*x*x-1.6*x))
#L-4:Distribution (Start Program)
#L10:PDF
#L11:CDF
#L12:Menu entry 2
#L13:Menu entry 3
#L14:Menu entra 4
#RI:Target for SOLVE
# ------------------------------------------------------------------------------

   000 {             } 
   001 {    42 21 14 } f LBL D
   002 {           1 } 1
   003 {    43 30  8 } g TEST x<y
   004 {    22 48  2 } GTO .2
   005 {           0 } 0
   006 {       44 25 } STO I
   007 {       43 33 } g R⬆
   008 {    32 48  0 } GSB .0
   009 {       43 32 } g RTN
   010 { 42 21 48  2 } f LBL .2
   011 {          33 } R⬇
   012 {           2 } 2
   013 {    43 30  8 } g TEST x<y
   014 {    22 48  3 } GTO .3
   015 {          33 } R⬇
   016 {          33 } R⬇
   017 {       44 25 } STO I
   018 {           2 } 2
   019 { 42 10 48  0 } f SOLVE .0
   020 {       43 32 } g RTN
   021 { 42 21 48  3 } f LBL .3
   022 {          33 } R⬇
   023 {           3 } 3
   024 {    43 30  8 } g TEST x<y
   025 {    22 48  4 } GTO .4
   026 {           0 } 0
   027 {       44 25 } STO I
   028 {       43 33 } g R⬆
   029 {    32 48  1 } GSB .1
   030 {       43 32 } g RTN
   031 { 42 21 48  4 } f LBL .4
   032 {          33 } R⬇
   033 {          33 } R⬇
   034 {       44 25 } STO I
   035 {           0 } 0
   036 { 42 10 48  1 } f SOLVE .1
   037 {       43 32 } g RTN
   038 { 42 21 48  0 } f LBL .0
   039 {       43 11 } g x²
   040 {           2 } 2
   041 {          16 } CHS
   042 {          10 } ÷
   043 {          12 } eˣ
   044 {       43 26 } g π
   045 {           2 } 2
   046 {          20 } ×
   047 {          11 } √x̅
   048 {          10 } ÷
   049 {       45 25 } RCL I
   050 {          30 } −
   051 {       43 32 } g RTN
   052 { 42 21 48  1 } f LBL .1
   053 {          36 } ENTER
   054 {          36 } ENTER
   055 {          36 } ENTER
   056 {          20 } ×
   057 {          20 } ×
   058 {          48 } .
   059 {           0 } 0
   060 {           7 } 7
   061 {          16 } CHS
   062 {          20 } ×
   063 {          34 } x↔y
   064 {           1 } 1
   065 {          48 } .
   066 {           6 } 6
   067 {          20 } ×
   068 {          30 } −
   069 {          12 } eˣ
   070 {           1 } 1
   071 {          40 } +
   072 {          15 } 1/x
   073 {       45 25 } RCL I
   074 {          30 } −
   075 {       43 32 } g RTN

# ------------------------------------------------------------------------------

Have fun!

Regards
deetee

[Albert Chan]

(11-24-2021 03:37 PM)deetee Wrote:  To calculate the CDF I used some approximating function rather than integrating the PDF,
which offers a accuracy of at least 2 digits (error < 1%) and is much faster.

More like 3 digits, and, formula very compact !

cdf2(z) := 1/(1+exp(-0.07*z^3-1.6*z))

We need an extra term to push it to 4-digits accuracy.

cdf3(z) := 1/(1+exp(0.0008*z^5-0.0743*z^3-1.595*z))

A Sigmoid Approximation of the Standard Normal Integral, by Gary R. Waissi and Donald F. Rossin
Note: link won't work as-is. Remove trailing characters after .pdf, and try again.

---

Approximation formula intentionally skipped fitting even powers of z-score.
This made formula matching identity: 1 - cdf(-z) = cdf(z)

1 - 1/(1+exp(-x)) = exp(-x)/(1+exp(-x)) = 1/(1+exp(x))

[deetee]

(11-25-2021 01:27 PM)Albert Chan Wrote:  More like 3 digits, and, formula very compact !

cdf2(z) := 1/(1+exp(-0.07*z^3-1.6*z))

We need an extra term to push it to 4-digits accuracy.

cdf3(z) := 1/(1+exp(0.0008*z^5-0.0743*z^3-1.595*z))

Thanks for sharing this information and link.

I developed my formula spending many hours with excel varying as few and simple parameters as possible.

After reading your document and doing some research I found a similar result from Bowling et al (2009) with a maximal error of 1.4x10-4:

CDF = 1/(1+exp(-0.07056*z^3-1.5976*z))

Despite there are many efforts from 1946 (Polya) to 2016 (Eidous and Al-Salman) it seems that no one has really solved this problem (simple formula, simple parameters, high accuracy) ...

Regards
deetee

[John Keith]

(11-28-2021 11:44 AM)deetee Wrote:  Despite there are many efforts from 1946 (Polya) to 2016 (Eidous and Al-Salman) it seems that no one has really solved this problem (simple formula, simple parameters, high accuracy) ...

Perhaps not, but numerical integration is plenty fast on modern hardware, if not on the 15C. I believe that the distribution functions on the 48 series and the Prime use numerical integration.

[Albert Chan]

(11-28-2021 11:44 AM)deetee Wrote:  I found a similar result from Bowling et al (2009) with a maximal error of 1.4x10-4:

CDF = 1/(1+exp(-0.07056*z^3-1.5976*z))

Nice ! And, we can easily invert for z, without SOLVE

ln(1/CDF-1) = -0.07056*z^3 - 1.5976*z       // already a depressed cubic

With cubic and linear term having same sign, use identity: sinh(3θ) = 4*sinh(θ)^3 + 3*sinh(θ)

→ z ≈ 5.494 * sinh(asinh(-0.3418*ln(1/CDF-1))/3)

---

We can use above approximation for the error function, and its inverse.

erfc(x) = 1 - erf(x) = cdf(-x*sqrt(2))*2 ≈ 2/(1+exp(0.2*x^3+2.259*x)

ierfc(x) = icdf(x/2)/-sqrt(2) ≈ 3.885 * sinh(asinh(0.3418*ln(2/x-1))/3)

--> ierf(x) ≈ 3.885 * sinh(asinh(0.6836*atanh(x))/3)

[deetee]

(11-28-2021 05:37 PM)Albert Chan Wrote:  → z ≈ 5.494 * sinh(asinh(-0.3418*ln(1/CDF-1))/3)

Hello Albert!

Thanks for your formula - but I have problems because CDF has to be between 0 and 1 - and CDF-1 is negative, where the ln ist not defined ... hmm (or am I wrong?).

But inspired by your "jaunty and cool " math-handling (and the TV-News from yesterday, where COVID-scientists worked with PDF and CDF functions) I developed some (new) formulae by myself (using hyperbolic functions):

CDF(x) ≈ (1 + tanh(7/200* x^3 + 4/5 * x)) / 2 ... Maximal error < 4.10-4

And the corresponding Inverse-CDF (one real root of this cubic equation):

x(CDF) = a * sinh(1/3 * asinh(b * atanh(1 - 2*CDF)))

where a = -2 * sqrt(160/7/3) ≈ -5.52
and b = 3/32 * sqrt(7*15/2) ≈ 0.679

Probably I don't change my original HP-15C-program because of many more program steps for a little bit faster calculation (no SOLVE) but for my standard linux console calculator (CLAC, FORTH-like, no SOLVE-function) it's worth to define new functions for PDF, CDF and their inverse.

Regards
deetee

[Albert Chan]

Hi, deetee

With 0 < CDF < 1, 1/CDF - 1 > 0

Also, cdf/icdf using tanh/atanh is equivalent to original exp/ln

(1 + tanh(x))/2 = (sinh(x)+cosh(x)) / (2*cosh(x)) = e^x / (e^x + e^-x) = 1 / (1 + e^(-2x))

CDF ≈ 1/(1+exp(-0.07056*z^3-1.5976*z)) = (1 + tanh(0.03528*z^3+0.7988*z))/2

Similarly, for the inverse, ln(x) = 2*atanh((x-1)/(x+1))

z ≈ 5.494 * sinh(asinh(-0.3418*ln(1/CDF-1))/3) = 5.494 * sinh(asinh(-0.6836*atanh(1-2*CDF))/3)

[deetee]

(11-30-2021 01:28 PM)Albert Chan Wrote:  With 0 < CDF < 1, 1/CDF - 1 > 0

Hi Albert!
Sorry, that was my fault. I read 1/CDF-1 as 1/(CDF-1)
Thanks for solving my comprehension problem.

Regards
deetee

[spartacus123111]

(11-28-2021 05:37 PM)Albert Chan Wrote:  

(11-28-2021 11:44 AM)deetee Wrote:  I found a similar result from Bowling et al (2009) with a maximal error of 1.4x10-4:

CDF = 1/(1+exp(-0.07056*z^3-1.5976*z))

Nice ! And, we can easily invert for z, without SOLVE

ln(1/CDF-1) = -0.07056*z^3 - 1.5976*z       // already a depressed cubic

With cubic and linear term having same sign, use identity: sinh(3θ) = 4*sinh(θ)^3 + 3*sinh(θ)

→ z ≈ 5.494 * sinh(asinh(-0.3418*ln(1/CDF-1))/3)

---

We can use above approximation for the error function, and its inverse.

erfc(x) = 1 - erf(x) = cdf(-x*sqrt(2))*2 ≈ 2/(1+exp(0.2*x^3+2.259*x)

ierfc(x) = icdf(x/2)/-sqrt(2) ≈ 3.885 * sinh(asinh(0.3418*ln(2/x-1))/3)

--> ierf(x) ≈ 3.885 * sinh(asinh(0.6836*atanh(x))/3)

Hi Albert,


Could yoou please explain how you arrived at:

z ≈ 5.494 * sinh(asinh(-0.3418*ln(1/CDF-1))/3)

I dont see how to get it even using the identity.

Thanks!

[AnnoyedOne]

FYI HP published a program for solving the Normal Distribution function on the HP-15C back in 1982. It uses numerical integration.

Page 60 of HP-15C Advanced Functions Handbook, August 1982

Page 51 of HP-15C Advanced Functions Handbook, HP 2011 reissue

Page 60 of HP-15C Collector's Edition Advanced Functions Handbook

A1

[Albert Chan]

(09-03-2024 03:45 PM)spartacus123111 Wrote:  

(11-28-2021 05:37 PM)Albert Chan Wrote:  Nice ! And, we can easily invert for z, without SOLVE

ln(1/CDF-1) = -0.07056*z^3 - 1.5976*z       // already a depressed cubic

With cubic and linear term having same sign, use identity: sinh(3θ) = 4*sinh(θ)^3 + 3*sinh(θ)

→ z ≈ 5.494 * sinh(asinh(-0.3418*ln(1/CDF-1))/3)

Could yoou please explain how you arrived at:

z ≈ 5.494 * sinh(asinh(-0.3418*ln(1/CDF-1))/3)

I dont see how to get it even using the identity.

Sorry, I missed your post ...

We just scale coefficients to match. First, scale z^3 coeff = 4

-56.6893*ln(1/CDF-1) = 4*z^3 + 90.5669*z

Next, scale z coeff to 3, ratio = 90.5669 / 3 = 30.1890 ≈ 5.494^2
Divide both side by 30.8190^(3/2) = 5.494^3

-0.3418*ln(1/CDF-1) = 4*(z/5.494)^3 + 3*(z/5.494)

This now match sinh(3θ) = 4*sinh(θ)^3 + 3*sinh(θ), solve for sinh(θ) = z/5.494

[Albert Chan]

(11-30-2021 01:28 PM)Albert Chan Wrote:  CDF ≈ 1/(1+exp(-0.07056*z^3-1.5976*z)) = (1 + tanh(0.03528*z^3+0.7988*z))/2

Similarly, for the inverse, ln(x) = 2*atanh((x-1)/(x+1))

z ≈ 5.494 * sinh(asinh(-0.3418*ln(1/CDF-1))/3) = 5.494 * sinh(asinh(-0.6836*atanh(1-2*CDF))/3)

If we go for central area, formula is simpler.

p = CDF(-z .. z) = 2 * CDF(0 .. z)

lua> p = fn'z: tanh(0.03528*z^3+0.7988*z)'
lua> z = fn'p: 5.494 * sinh(asinh(0.6836*atanh(p))/3)'

lua> p(1), p(2), p(3)
0.6826606943611349      0.9544778797318427      0.9975360894437647