Simulate an M/M/1 queue and calculate moments of the distribution of total running times of k jobs encountered by an arriving job; see http://math.stackexchange.com/questions/218234

Question218234.java

import java.util.Random;

public class Question218234 {
    final static Random random = new Random ();

    static double exponential (double mu) {
	return -Math.log (random.nextDouble ()) / mu;
    }

    final static double lambda = 2;
    final static double mu = 3;

    final static int k = 2;

    public static void main (String [] args) {
	int queue = 0;

	double time = 0;
	double nextArrival = exponential (lambda);
	double nextDeparture = 0;
	double currentDuration = 0;

	long total = 0;
	double sum = 0;
	double sum2 = 0;

	while (time < 100000000) {
	    if (queue == 0 || nextArrival < nextDeparture) {
		if (queue == k) {
		    double duration = currentDuration;
		    for (int i = 1;i < k;i++)
			duration += exponential (mu);
		    total++;
		    sum += duration;
		    sum2 += duration * duration;
		} 
		time += nextArrival;
		if (queue != 0)
		    nextDeparture -= nextArrival;
		nextArrival = exponential (lambda);
		queue++;
	    }
	    else {
		time += nextDeparture;
		nextArrival -= nextDeparture;
		nextDeparture = currentDuration = exponential (mu);
		queue--;
	    }
	}

	System.out.println (sum / total + " / " + sum2 / total);
    }
}