HP Forums
(15C)(DM15) - PDF/CDF and Inverse of a Normal Distribution - Printable Version

+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Software Libraries (/forum-10.html)
+--- Forum: General Software Library (/forum-13.html)
+--- Thread: (15C)(DM15) - PDF/CDF and Inverse of a Normal Distribution (/thread-17744.html)



(15C)(DM15) - PDF/CDF and Inverse of a Normal Distribution - deetee - 11-24-2021 03:37 PM

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


RE: (15C)(DM15) - PDF/CDF and Inverse of a Normal Distribution - Albert Chan - 11-25-2021 01:27 PM

(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) Smile

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


RE: (15C)(DM15) - PDF/CDF and Inverse of a Normal Distribution - deetee - 11-28-2021 11:44 AM

(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


RE: (15C)(DM15) - PDF/CDF and Inverse of a Normal Distribution - John Keith - 11-28-2021 02:02 PM

(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.


RE: (15C)(DM15) - PDF/CDF and Inverse of a Normal Distribution - Albert Chan - 11-28-2021 05:37 PM

(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)


RE: (15C)(DM15) - PDF/CDF and Inverse of a Normal Distribution - deetee - 11-30-2021 11:07 AM

(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


RE: (15C)(DM15) - PDF/CDF and Inverse of a Normal Distribution - Albert Chan - 11-30-2021 01:28 PM

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)


RE: (15C)(DM15) - PDF/CDF and Inverse of a Normal Distribution - deetee - 11-30-2021 02:46 PM

(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