Origins of HP41 numerical routines ot compute statistical distributions

[Les Wright]

This is more of a theoretical question.

In the Stat Pac module there is a routine, SigmaNormd, that computes a number of features relevant to the standard normal distribution--i.e., actual value of the density function for a given z, upper-tail prob of given z, z of a given upper tail probability.

On inspecting the code, it looks like the probability density function is approximated algebraically in some way--for example, it looks like there are a number of numeric constants stored in the first few registers on initialization of the routine, and these get recalled later as the relevant estimates are computed. In other words, it doesn't look as though the program works directly with the actual density function (i.e., (1/sqrt(2*Pi))*exp(-z^2/2)).

Does anyone know the origins of this routine? Is it from Numerical Recipes? Is it based on a Taylor series approximation? (As for the latter, I have found on my HP48 that one has to go to a pretty high order in order for the polynomial approximation to maintain accuracy in the upper and lower tails. In emulation this is fast enough, but on the actual calculator its as slow as molasses.)

This is more of a point of curiosity than of great practical relevance, but if anyone knows what the actual formulae used are, and their source, I would be much obliged.

Les

[Les Wright]

Just posting to correct typo in subject line.

[Joao]

From what you say, it looks like a polynomial approximation of the type described by

C. Hastings, Jr. Approximations for Digital Computers, Princeton University Press, 1955.

Can you tell us the values of the constants?

Please let us know if you find more about this.

Joao

Edited: 8 May 2006, 5:07 p.m.

[Les Wright]

I won't bore you with all the digits.

From what I can see, the routine has some similarities to a nested routine in Numerical Recipes to compute the complementary error function, erfc, which of course is intimately related to the cumulative normal distribution. The coefficients used there are not the same, but that doesn't matter to me--the bottom line is that I have discerned that there is something a bit more complex going on than a simple Taylor/Maclaurin series expansion.

Ah, I wish I knew more about approximation theory--Chebyshev polynomials, Fourier series, etc.

Les

[Les Bell]

Les Wright wrote:

Quote:

---

Ah, I wish I knew more about approximation theory--Chebyshev polynomials, Fourier series, etc.

---

Les, take a look at http://www.google.com.au/search?q=Scientific+Analysis+on+the+Pocket+Calculator+Jon+Smith&start=0&start=0&ie=utf-8&oe=utf-8&client=firefox-a&rls=org.mozilla:en-US:official . I notice one went on eBay just recently for only $5.

Best,

--- Les Bell

[http://www.lesbell.com.au]

[Trent Moseley]

I can't help but to respond to your last posting on Jon M. Smith's book "Scientific Analysis on the Pocket Calculator". It rang a bell (no pun to you). I hadn't used it in years, but I found it in a dusty corner, 1977 2nd ed. And gee it's only worth five bucks.

tm

[Les Bell]

Yes - I have a copy on the bookshelf, myself. Bought it in 1976, just after I got my '45. It was very useful!

Best,

--- Les Bell
[http://www.lesbell.com.au]

[Les Wright]

Actually, the code does use the exact function for determining the ordinate of the actual density function, but seems to use an interesting approximation when dealing with the cumulative distribution, since it is related to the error function erf(x) which does not have an exact indefinite integral.

I totally forgot that I have both Numerical Recipes in Pascal (the book) and Numerical Recipes in C (the code files), both purchased some years ago. The following routine estimates complementary error function erfc(x), for nonnegative x. For negative x, compute erfc(|x|) and subtract from 2. To get erf(x) subtract erfc(x), whether x is positive or negative, from unity. The relationship between erf(x) and the cumulative normal dist. has to do with scaling factors placed upon the argument and the overall function itself.

Here is that code. For those who don't recognize the syntax, it is cut and paste from Maple:

z:=abs(x);
t:=1/(1+0.5*z);
ans:=t*exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+t*(0.09678418+t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+t*(-0.82215223+t*0.17087277)))))))));

The constants and the routine in the code I was wondering about doesn't look quite like this, BUT I can discern that something like this is going on.

I will gladly look into that classic book recommended in this thread.

Les

[Mike (Stgt)]

Your "density function" (1/sqrt(2*Pi))*exp(-z^2/2) may be only be integrated anlytically by squaring the function, integrating and takeing the square root of the result. Integrating numerically may be converging too slow. More in here.

Ciao.....Mike

[Namir]

The equations to calculate the popular probability distribution functions and their inverse are found in the Hadbook Of Mathematical Functions.

In my book "Mathematical Algorithms in Visual Basic for Scientists & Engineers" I show you libraries that calculate these functions and their inverses. You can still find a copy of this book in Amazon. Search for "Used Book" after you select the Books category, since the book is about 10 years old. In my book I show you the pseudo-code AND the actual code in Visual Basic". I also wrote a C+++ version of the book.

Namir Shammas

[Les Wright]

Though I understand I cannot identify the source, I did take a look at the Business Stats solution book on a certain website. It uses a similar routine (at least the code looks similar) to that in the Finance Pac, and quotes the Handbook of Mathematical Functions as the source. So too does Numerical Recipes.

The Handbook is in a Dover reprint and is readily available in Canada for under $40, so maybe I will get it new or used. I will all keep a look out for the Bell book (one eBay seller offers it for $2 but doesn't ship outside of the US) and Namir's book.

Les

[Namir]

Les,

Check a orevious (or archived) post. Someone on this site has posted a link to the PDF version of the Handbook of Mathematical Functions.

Namir