public static int poisson(double lambda) {
        if (lambda <= 0.0)
            throw new IllegalArgumentException("Parameter lambda must be positive");
        // using algorithm given by Knuth
        // see http://en.wikipedia.org/wiki/Poisson_distribution
        int k = 0;
        double p = 1.0;
        double L = Math.exp(-lambda);
        do {
            k++;
            p *= uniform();
        } while (p >= L);
        return k-1;
    }