I would like to create a function in Java that generates Poisson arrivals given the mean arrival rate (lambda) and the mean service rate (mu).
In my example I have: 2,2 requests/day, in other words 2,2 arrivals/day and a mean service time of 108 hours. Considering that my program starts at t=0 minutes, I would like to create a function that returns arrivals[], which will contain, t1, t2, and a possible t3. T1,t2 and t3 are the instants (in minutes) during the day where this arrivals occur. I have the following restrictions:
t1 < t2 < t3 < 1440 minutes (24 hours*60 minutes/hour)
t2-t1 > 108 minutes
t3-t2 > 108 minutes
t3+ 108 minutes < 1440 minutes
Can someone please help me?
Thank you,
Ana
You can use this algorithm proposed by D. Knuth:
private static int getPoissonRandom(double mean) {
Random r = new Random();
double L = Math.exp(-mean);
int k = 0;
double p = 1.0;
do {
p = p * r.nextDouble();
k++;
} while (p > L);
return k - 1;
}
To understand how this works note that after k iterations the loop condition becomes
p1 * p2 * ... * pk > L
which is equivalent to
-ln(p1)/mean -ln(p2)/mean ... -ln(pk)/mean > 1
Note that if p is uniformly distributed then -ln(p)/mean has exponential distribution with a given mean. A random variable with Poisson distribution is equal to the number of times a given event occurs within a fixed interval when the lengths of the intervals between events are independent random variables with exponential distribution. Since we use the mean of the Poisson distribution as the mean of the exponential distribution for the intervals between events, the fixed internal in which we count occurrences is unit length. Thus, the loop condition sums up the lengths of the intervals between events and checks whether we have gone beyond the unit interval. If we have gone beyond the unit interval while counting kth event, then k-1 events occurred within the interval and so we return k-1.
I found this solution, using inverse transform sampling:
http://preshing.com/20111007/how-to-generate-random-timings-for-a-poisson-process
It doesn't use a rejection sampling approach, is efficient and accurate.
It uses the fact that the distribution of the time between events is an exponential distribution, with parameter lambda which is the arrival rate.
The Exponential distribution is lambda exp(-lambda x).
In order to sample values from that distribution and avoid rejection sampling, it is easier to use its Cumulative Distribution Function (CDF): 1 - exp(-lambda x).
The CDF is a function that starts at 0.0 and grows to 1.0 with larget parameters.
Intuitively, the probability that you will get an event increases as time passes.
To sample the CDF, and again avoid rejection sampling, it is easier to pick a uniform random number U between [0,1] and to plug that value in the inverse function of the CDF, which gives: nextEvent = - Ln(U)/lambda.
Because Ln(0) is undefined, and most random number generators include 0.0 and exclude 1.0, it is safer to use:
nextEvent = -Ln(1.0-U)/lambda.
If your code uses a millisecond based time/sleep function you can use:
double rate = 20.0/1000.0; // 20 per second on average
sleep(floor(-1.0 * log(1.0 - rand()*1.0/RAND_MAX) / rate) );
you can add this to build.gradle
implementation 'org.kie.modules:org-apache-commons-math:6.5.0.Final'
and use class PoissonDistribution more detail
Here's some simplified code to generate a poisson number with a given mean:
private static int poisson(double mean) {
int r = 0;
double a = random.nextDouble();
double p = Math.exp(-mean);
while (a > p) {
r++;
a = a - p;
p = p * mean / r;
}
return r;
}
You should be able to use this or something similar to generate your arrival numbers for each time period: the input should be the mean number of arrivals you expect in the period.
Note that if your mean is very large (say 500+) you will want to approximate the number of arrivals with a normal distribution instead. This is more efficient, plus it avoid numerical overflow issues inherent in the code above (at some point Math.exp(-mean) gets rounded to zero......ooops!)
