+- 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)
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)
--> ierf(x) ≈ 3.885 * sinh(asinh(0.6836*atanh(x))/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
RE: (15C)(DM15) - PDF/CDF and Inverse of a Normal Distribution - spartacus123111 - 09-03-2024 03:45 PM
(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!
RE: (15C)(DM15) - PDF/CDF and Inverse of a Normal Distribution - AnnoyedOne - 09-04-2024 02:18 PM
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
RE: (15C)(DM15) - PDF/CDF and Inverse of a Normal Distribution - Albert Chan - 10-28-2024 11:22 PM
(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
RE: (15C)(DM15) - PDF/CDF and Inverse of a Normal Distribution - Albert Chan - 10-30-2024 03:49 PM
(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
![[Image: empirical-rule-normal-distribution-graph.png]](https://www.inchcalculator.com/wp-content/uploads/2021/11/empirical-rule-normal-distribution-graph.png)