# HP Forums

You're currently viewing a stripped down version of our content. View the full version with proper formatting.
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
#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
(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))
(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
(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.
(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)
(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
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)
(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
Reference URL's
• HP Forums: https://www.hpmuseum.org/forum/index.php
• :