The Museum of HP Calculators

HP Forum Archive 21

[ Return to Index | Top of Index ]

One for Mr NAMIR
Message #1 Posted by PGILLET on 22 Mar 2013, 4:51 p.m.

Here is a short HP15C (DM15CC) program for normal distribution
(See below for source)

********** Phi(X) **********
LBL A
STO 2
STO 3
X^2
STO 4
1
STO 5
0
STO 6

LBL 1
2
STO+ 5
RCL 4
RCL/ 5
STO* 3
RCL 2
STO 6
RCL+ 3
STO 2
RCL- 6
ABS
1E-6
-
TEST 3
GTO 1
RCL 4
2
/
CHS
9.189385332E-1
-
E^X
RCL* 2
2
1/X
+
RTN

********** Phi-1(X) **********
LBL B
STO 1
SOLVE 2
RTN

LBL 2
GSB A
RCL- 1
RTN


Phi(X) based on this powerful algorithm:
Marsaglia in 2004:
http://www.jstatsoft.org/v11/a05/paper
http://stat-athens.aueb.gr/~jpan/papers/Panaretos-JIRSS2008%2857-72%29ft.pdf

public static double normalCDF(double x) {
double s = x, t = 0., b = x, q = x*x, i = 1.;
while(s != t) {
i += 2.;
b = b * (q/i);
t = s;
s = t + b;
}
return .5 + s * Math.exp(-.5 * q - .91893853320467274178);
}

.91893853320467274178 is in fact LN(2*PI)/2

Register usage in Phi(X) routine:
R2: s
R3: b
R4: q
R5: i
R6: t

Instead of testing s=t for ending Phi(x), I preferred testing abs(s-t)<10^-6, in order to tune the accuracy and response time:
Phi(1): 0.8413447388 (1 sec) (OK 7 significants digits)
Phi(5): 0.9999997153 (3.5 sec)(OK 8 significant digits)
Phi-1(.995): 2.575829318 (25 sec) (OK 8 significant digits)
Phi-1(.99995): 3.890588721 (49 sec) (OK 5 significant digits)

      
Re: One for Mr NAMIR
Message #2 Posted by Namir on 23 Mar 2013, 5:13 a.m.,
in response to message #1 by PGILLET

Dear PGillet.

I am 110% sincere when I say that I read you message with amazement and a great sense of thank you!. Recently, I have been looking at a statistical distribution handbook that includes pseudo-code to generate random numbers that are uniform, normal, chi-distributed, Student-t distributed, and F distributed. I was looking at these routines as good short programs to use in playing with my (now) vast collection of calculators and BASIC pocket PCs. What you provided me is, as we say, right on the money!

I am praising you sense of psychic connection of your unconscious mind in detecting that I needed a routine like the one you shared!!!

So, thank you, thank you, thank you!! I have already saved the information in your message.

Namir

PS: If you have the time and the inkling, can you locate code for random numbers for Student, Chi-Square, and F distributions.

Edited: 23 Mar 2013, 9:24 a.m.

            
Re: One for Mr NAMIR
Message #3 Posted by Dale Reed on 23 Mar 2013, 5:16 p.m.,
in response to message #2 by Namir

Namir and all,

I'm sure all you mathematician types are aware (but I wasn't until not long ago, not being a math major) of the Box-Muller Transform.

Here's the Wikipedia article about it. I use this method with industrial programmable controllers to simulate normally-distributed noise (as would be seen on a sensor input) to test analog input processing, as well as PID and other advanced process control methods.

It's super easy to use in such an environment. Two random number calls and two calculation blocks. (The controller has built-in sin, cos, tan, sqrt, etc., for floating-point.) I ran some logic for a few weeks (several tens of millions of samples) that captured how many noise points were between 0.0 and 0.1, between 0.1 and 0.2, etc., from <-3.0 to >+3.0 and the curve was right on, except for a tiny "dent" around 0.1<= noise <0.2.

I also kept a circular buffer of the last 100 points and ran mean and std. The mean stays very close to zero, and standard deviation stays very close to 1.0.

I attribute any imperfections to the random number generator (something I made that's not guaranteed to be perfectly uniformly distributed) and to the use of single-precision float (32-bit).

Dale

                  
Re: One for Mr NAMIR
Message #4 Posted by Dieter on 24 Mar 2013, 10:15 a.m.,
in response to message #3 by Dale Reed

Dale,

Quote:
I'm sure all you mathematician types are aware (but I wasn't until not long ago, not being a math major) of the Box-Muller Transform.
While I am certainly not a mathematician type, I heard about the Box-Muller method (without knowing its name) more than 30 years ago when I read the manual of TI's standard library module (Prgm 15 provides a random number generator). At that time I did not know much about statistical distributions, and I wondered why the given formula contained a trig function - which in my litte world appeared only in geometry. #-)

There are other (also faster) methods for generating normally-distributed random numbers. A very simple one that is often "good enough" simply adds twelve evenly distributed random numbers and decreases this sum by six. No trigs, no logs, no roots required. ;-)

Dieter

      
Re: One for Mr NAMIR
Message #5 Posted by Dieter on 23 Mar 2013, 9:59 a.m.,
in response to message #1 by PGILLET

The method used here seems to be the well-known series expansion for the Normal integral (e.g. as in Abramovitz & Stegun 26.2.11). It works fine for the central region, but less so far out in the tails. That's why the Marsaglia paper suggests a different approach there. Personally, I prefer a combination of this series expansion up to |z| = 2,326 (where Phi drops below 0,01) and a continued fraction approach (e.g. A+S 26.2.14) with a pre-calculated number of iterations. This is also the way the 34s originally used to evaluate the Normal integral. It is both fast and accurate.

Regarding the inverse (quantile) function: if you want to use the solver I would suggest to provide two decent initial guesses. This is easily accomplished: let's assume p =< 0.5 (the other half is handled via quantile(1-z) = -quantile(z)), then an upper bound for the quantile is sqrt(-2 ln p). A lower bound then simply can be defined by the upper bound minus sqrt(ln 4). This should substantially speed up the Solver.

I would also like to suggest a look at the HP15C Advanced Functions Handbook, p. 51 ff. It shows another approach for the Normal integral, this time evaluated with the 15C's Integrate function combined with a clever variable substitution.

By the way - I'm a big fan of rational approximations. If the goal is merely six significant digits over the central range (say, 1E-10 < p < 1-1E-10) designing such an approximation is quite trivial. It just requires a number of pre-stored constants. Results both for the CDF and the quantile should then be returned almost instantly. ;-)

The Normal distribution (both the CDF and the quantile function) has been discussed here several times with regard to their implementation for the 34s project. You will find several threads on this topic here.

Dieter

[ Return to Index | Top of Index ]

Go back to the main exhibit hall