Is there a simple and efficient way of implementing a function for the standard cumulative (standard) normal distribution with the Solve utility?
I'm aware of the Abramowitz polynomial function but maybe there is a more efficient way like a series expansion I could use?
Grateful for any advice. I'm using an HP-27S.
Rgds, Anders
I don't have Abramowitz and Stegun in front of me but I think that a rational function approximation is the way to go. I don't know how many program steps you have in an HP27S, but if the Solver works as I think it should you could use the rational function approximation for the cumulative normal distribution - percentile = 0 as your main equation, and for the inverse distribution (finding a particular z value for a given percentile or alpha) use the solver--assuming that that is what the Solver does. I am not familiar with the calculator.
The HP41CV routine NORMAL (http://www.hpmuseum.org/software/41td/normdst.htm) does it all but requires over 160 steps. A lot of this is devoted to storing the various coefficients in registers and comprises estimations for both of the cumulative normal and inverse normal distributions, so with a Solver utility you could pare out a lot of these steps.
I have played around with straight forward integrations of Maclaurin series expansions of the normal curve function on the HP48G and they take many terms to converge, particularly in the tails. Simpson's integration of the curve (1/sqrt(2*Pi) * exp(-t^2/2) using INTG from the HP41 MathPac module works too but is slow if you take a reasonable number or sampling points.
This all said,I could be talking out of my ear here--I know nothing about the strengths and limitations of the HP27S.
Any other thoughts out there?
Les
p.s. I had a glancingly related question regarding normal distribution calculations a few weeks ago. On the groups suggestion I got the Dover reprint of Abramowitz and Stegun, and am glad I did. An oldie but a goodie. What a treasure trove!
I now have Abramowitz and Stegun in front of me, and on closer inspection I can recommend the approximation of the inverse distribution, 26.2.23, with only modest enthusiasm. Absolute error quoted as less than 4.5e-4, which means the approximation is really only okay to 3 or at most 4 decimal places. I have run the code I referenced in the last post under P41CX, V41 and Free42, and, sure enough, the z-score associated with key alphas such as .01, .05, .001, etc matches more precise (? accurate) comparators such as UTPN in the HP48 series or the relevant computations in Maple up to about the fourth decimal place at best.
On the other hand, 26.2.17 in Abramowitz and Stegun claims absolute error less than 7.5e-8. If you could implement this approximation as the backbone of your equation in a Solver, then solve for the desired z for a given percentile or alpha (i.e., 1-percentile), you may get more decimal places, hopefully correct ones.
What exactly are you hoping for? Tell me more about program memory of the HP27.
Les
The HP27S has the same type of solver as the HP17BII, HP18C, HP19B, etc.
It's a standard formula solver. IIRC, it does not have let and get, however.
The description at http://www.hpmuseum.org/hp27s.htm helps me out a bit. Looks like the calc has a solver similar to the one I am most familiar with in HP48G, albeit not as elaborate and permitting only algebraic expressions.
That clarified, I suggested to the original poster that he try approximation 26.2.17. If A(z) is the the approximation, and a is the alpha, then setting up the solver with A(z) - 1 + a = 0 is a good starting point. I should play around with the approximation in my hp48g and see where I end up. The only problem I see is that, with required substitutions, formula 26.2.17 can get pretty messy to key into the calculator.
Les
You can store the constants in a list called "NORM" and then use this:
{utpn=SIGMA(I:1:5:1:ITEM(NORM:I)*SPPV(23.16419*x:I))
*EXP(-x*x/2)}
or try this
{ndist=.5+sigma(n:0:30:(2*x^2)^n*FACT(n)/FACT(2*n+1))*x*EXP(-
x^2/2)/SQRT(2*PI)}
many others, and NORM (5 numbers adding to 1) in wronski.zip - attached to a DataFile article from a year ago - see www.hpcc.org
Cheers,
Tony
Neat!
In the first case you make clever use of the SPPV function, which I gather is a financial function found on this calculator that just happens to generate the desired transformation of the independent variable.
In the second case you use a series expansion out to 30 terms that capitalizes on the relationship between the cumulative normal distribution and the incomplete gamma function, which has a series expansion that converges rapidly for smallish arguments. On my hp48 this gives excellent accuracy when compared with the calculator's intrinsic UTPN function. The drawback is that when I use it in the equation solver to find the Z associated with a given probability, the sucker chugs along for ever as the calculator calculates the series expansion for every trial value of Z. Is the equation solver in the HP27S pretty swift? You could speed things up by taking the expansion out to fewer terms, but you lose decimal accuracy, especially as Z gets bigger and you get more in tail of the bell curve. (Please not that I use Z here to refer to a standard normal variate, not a complex number.)
It will be interesting to see what the original poster decides to do. Anders, can you let us know where your exploration takes you?
Les
Hi again,
many thanks for your efforts in trying to help me out on this! I tried the
{ndist=.5+sigma(n:0:30:(2*x^2)^n*FACT(n)/FACT(2*n+1))*x*EXP(- x^2/2)/SQRT(2*PI)}
suggested by Tony and it works well for my needs. In fact, since I only need to evaluate NDIST(x) for -2<x<2, I can decrease the summation from 30 to 10 steps and still get 4-5 decimals accuracy, and a bit faster calculation.
The only drawback of this method is that it doesn't seem to work for x=0?
I don't have the actual Abramowitz book (should probably get it), but I guess 26.2.17 equation is the "standard" polynominal often referenced? Haven't been able to test it, though.
Once more, thanks
//Anders
Dead right Les - the 27S has SPPV (I checked on Craig Finseth's page before posting<G>), and the 27S solver wouldn't be swift. The 12c has NPV which is programmable to boot and with the constants stored as cashflows one has a very very economical UTPN - it really surprised me.
Cheers,
Tony
Hi Anders - what do you get for x=0? It should give ndist=.5?
Or is it a problem when re-solving for x?
Yup the Abramowitz 26.2.17 is often used. The first equation I gave does that - you just need the 5 constants in the list NORM.
But your 10 term ndist sounds pretty nifty! Thanks for the feedback. Glad you got it going.
Cheers,
Tony
I don't know enough about this (yet) to know what the "standard" polynomial is, but 26.2.17 is:
P(x) is approximated by 1 - Z(x)(b1*t + b2*t^2 ... b5*t^5),
where t = 1/(1 + px), P(x) is the cumulative normal distribution percentile (1 - upper tail alpha), Z(x) is the standard gaussian density function exp(-x^2/2)/SQRT(2*Pi), and p and the 5 b's are certain coefficients that I am too lazy to type now, since each includes 7 to 10 digits ;).
Abramowitz and Stegun gives three other examples that are simpler but with more reported error.
The series expansion should give .5 for x = 0. Did you enter the code correctly? If you plug in x = 0, all the terms except the leading .5 constant should vanish. Works fine on my HP48G. If you want to use values of x less than 0, you have to use the absolute value and subtract the result from 1--so NDIST(1) = about 0.84, but NDIST (-1) = about 0.16.
I have a question for Tony: did he leave out the 1/SQRT(2*Pi) accidentally in his code, or is it presumed factored into the above-mentioned b's?
Les
I should mention that I use Power48 and Emu48 more than my actual HP48G.
With "true calculator speed" turned OFF, in these emulators the solver speed is really quite impressive, even with 50 to 100 terms in the series expansion.
Another reason why I love emulators/simulators (see my recent post).
Les
Hi again Tony,
I just get 'Solution not found' when I set x=0. I don't understand why, as you say the second righthand term should be zero and make ndist=0.5
Anyway, it can easily be handled by an if statement.
An additional question: I guess there is no way of calling one user-defined function from within another function? Like if I want to use the ndist-function in two other functions?
Cheers, Anders
Les,
Yes, I had to factor the 1/SQRT(2*Pi) into the b's, as the 12C doesn't have PI :-) That's what makes the b's I use add to 1 exactly. I don't have the 5 b's to hand - they are in the Datafile article though.
Cheers,
Tony
Edited: 25 May 2006, 7:18 p.m.
Abramowitz and Stegun is online (see most recent thread started by Mr. Abillo). Here is the page with 26.2.17, and other goodies:
http://www.math.sfu.ca/~cbm/aands/page_932.htm