One for Mr NAMIR



#6

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)


#7

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.


#8

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


#9

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

#10

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


Possibly Related Threads...
Thread Author Replies Views Last Post
  41-MCODE: Dr. Jekyll & Mr. Hyde Ángel Martin 9 395 07-09-2012, 09:41 AM
Last Post: Monte Dalrymple
  For Mr. Diego Diaz Ignazio Cara (Italy) 1 111 04-30-2008, 03:29 PM
Last Post: Diego Diaz
  Mr.Hurd viz-a-viz the forthcoming HP-35s Trent Moseley 2 150 05-30-2007, 04:39 PM
Last Post: Eric Smith
  "Simple" Math Question (to keep Namir busy) Chuck 7 303 05-03-2007, 06:43 PM
Last Post: Peter A. Gebhardt
  Looking for Mr. David Smith. Was: Expiring EPROM's inHP's? Wolfgang Miller 3 209 09-09-2006, 07:45 PM
Last Post: David Smith
  To Mr. Eric Smith - about HP 9872A plotter Artur-Brazil 1 156 07-31-2006, 02:33 PM
Last Post: Eric Smith
  News from Mr Packard jose goncalves (Brasil) 6 274 07-16-2003, 07:30 AM
Last Post: Ron Ross
  Mr. Hicks how about a screen saver Norm 18 622 06-12-2003, 06:01 AM
Last Post: Meindert Kuipers
  Hello Mr. Sus Mike (Stgt) 0 96 05-21-2003, 03:48 AM
Last Post: Mike (Stgt)

Forum Jump: