CODEWITHSUNDEEP
javacurvepoisson
Dixon Steel
Dixon Steel

Reputation: 1061

Poisson Distribution in Java (correctness?)

I have to generate data for a Poisson distribution. My range is n = 1000 up to 100K. Where n is the number of data elements; k varies from 1 to n. It says to use lambda as n/2

I have never taken stats and have no idea how to get the correct curve here. I can feed it lambda as n/2, but do I vary K from 0-n? I tried this (passing k in as a parameter) and when I graphed the data it ramped up, not a fish tail. What am I doing wrong, or am I doing it correctly?

Thanks

I have this code in java from Knuth.

static double poissonRandomNumber(int lambda) {
    double L = Math.exp(-lambda);
    int k = 0;
    double p = 1;
    do {
        k = k + 1;
        double u = Math.random();
        p = p * u;
    } while (p > L);
    return k - 1;
}

Upvotes: 2

Views: 1528

Answers (1)

Jay Elston
Jay Elston

Reputation: 2078

One of the problems you are running into is a basic limitation of how computers represent and perform calculations with floating point numbers.

A real number is represented on a computer in a form similar to scientific notation:

Significant digits × base^exponent

For double precision numbers, there are 11 bits used for the exponent and 52 for the "significant digits" portion. Because floating point numbers are normalized, the first positive floating point number > 0.0 has a value of about 10^-320 (this is defined as Double.MIN_VALUE in Java). See IEEE Standard 754 Floating Point Numbers for a good writeup on this.

Consider the line of code: 

double L = Math.exp(-lambda);

With a lambda of 1000, e^-1000 (which is about 10^-435) is less than Double.MIN_VALUE, and there is no way the computer can represent e^-1000 any differently than it can represent e^-100000

You can solve this problem by noticing that lambda is an "arrival rate", and you can calculate random samples for shorter intervals and sum them. That is

x = p(L);

can be computed as

x = p(L/2) + p(L/2);

and larger numbers can be approximated:

x = 100 * p(L/100);

The Wikipedia article has on the Poisson distribution has some good pointers to ways to compute Poisson distributions for large values of lambda.

Upvotes: 1

Related Questions