PoissonSlow.java

/*
  Copyright � 1999 CERN - European Organization for Nuclear Research.
  Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose 
  is hereby granted without fee, provided that the above copyright notice appear in all copies and 
  that both that copyright notice and this permission notice appear in supporting documentation. 
  CERN makes no representations about the suitability of this software for any purpose. 
  It is provided "as is" without expressed or implied warranty.
*/
package sim.util.distribution;
import ec.util.MersenneTwisterFast;

/**
 * Poisson distribution; See the <A HREF="http://www.cern.ch/RD11/rkb/AN16pp/node208.html#SECTION0002080000000000000000"> math definition</A>
 * and <A HREF="http://www.statsoft.com/textbook/glosp.html#Poisson Distribution"> animated definition</A>.
 * <p>
 * <tt>p(k) = (mean^k / k!) * exp(-mean)</tt> for <tt>k >= 0</tt>.
 * <p>
 * Valid parameter ranges: <tt>mean > 0</tt>.
 * Note: if <tt>mean <= 0.0</tt> then always returns zero.
 * <p>
 * Instance methods operate on a user supplied uniform random number generator; they are unsynchronized.
 * <dt>
 * Static methods operate on a default uniform random number generator; they are synchronized.
 * <p>
 * <b>Implementation:</b> 
 * This is a port of <A HREF="http://wwwinfo.cern.ch/asd/lhc++/clhep/manual/RefGuide/Random/RandPoisson.html">RandPoisson</A> used in <A HREF="http://wwwinfo.cern.ch/asd/lhc++/clhep">CLHEP 1.4.0</A> (C++).
 * CLHEP's implementation, in turn, is based upon "W.H.Press et al., Numerical Recipes in C, Second Edition".
 *
 * @author wolfgang.hoschek@cern.ch
 * @version 1.0, 09/24/99
 */
public class PoissonSlow extends AbstractDiscreteDistribution {
    private static final long serialVersionUID = 1;

    protected double mean;

    // precomputed and cached values (for performance only)
    protected double cached_sq;
    protected double cached_alxm;
    protected double cached_g;

    protected static final double MEAN_MAX = Integer.MAX_VALUE; // for all means larger than that, we don't try to compute a poisson deviation, but return the mean.
    protected static final double SWITCH_MEAN = 12.0; // switch from method A to method B
        
    protected static final double[] cof = { // for method logGamma() 
        76.18009172947146,-86.50532032941677,
        24.01409824083091, -1.231739572450155,
        0.1208650973866179e-2, -0.5395239384953e-5};

/**
 * Constructs a poisson distribution.
 * Example: mean=1.0.
 */
    public PoissonSlow(double mean, MersenneTwisterFast randomGenerator) {
        setRandomGenerator(randomGenerator);
        setMean(mean);
        }
/**
 * Returns the value ln(Gamma(xx) for xx > 0.  Full accuracy is obtained for 
 * xx > 1. For 0 < xx < 1. the reflection formula (6.1.4) can be used first.
 * (Adapted from Numerical Recipes in C)
 */
    public static double logGamma(double xx) {
        double x = xx - 1.0;
        double tmp = x + 5.5;
        tmp -= (x + 0.5) * Math.log(tmp);
        double ser = 1.000000000190015;

        double[] coeff = cof;
        for (int j = 0; j <= 5; j++ ) {
            x++;
            ser += coeff[j]/x;
            }
        return -tmp + Math.log(2.5066282746310005*ser);
        }
/**
 * Returns a random number from the distribution.
 */
    public int nextInt() {
        return nextInt(this.mean);
        }
/**
 * Returns a random number from the distribution; bypasses the internal state.
 */
    private int nextInt(double theMean) {
        /* 
         * Adapted from "Numerical Recipes in C".
         */
        double xm = theMean;
        double g = this.cached_g;

        if (xm == -1.0 ) return 0; // not defined
        if (xm < SWITCH_MEAN ) {
            int poisson = -1;
            double product = 1;
            do {
                poisson++;
                product *= randomGenerator.nextDouble();
                } while ( product >= g );
            // bug in CLHEP 1.4.0: was "} while ( product > g );"
            return poisson;
            }
        else if (xm < MEAN_MAX ) {
            double t;
            double em;
            double sq = this.cached_sq;
            double alxm = this.cached_alxm;

            MersenneTwisterFast rand = this.randomGenerator;
            do { 
                double y;
                do {
                    y = Math.tan(Math.PI*rand.nextDouble());
                    em = sq*y + xm;
                    } while (em < 0.0);
                em = (double) (int)(em); // faster than em = Math.floor(em); (em>=0.0)
                t = 0.9*(1.0 + y*y)* Math.exp(em*alxm - logGamma(em + 1.0) - g);
                } while (rand.nextDouble() > t);
            return (int) em;
            }
        else { // mean is too large
            return (int) xm;
            }
        }
/**
 * Returns a random number from the distribution.
 */
    protected int nextIntSlow() {
        final double bound = Math.exp( - mean);
        int count = 0;
        double product;
        for (product = 1.0; product >= bound && product > 0.0; count++) {
            product *= randomGenerator.nextDouble();
            }
        if (product<=0.0 && bound>0.0) return (int) Math.round(mean); // detected endless loop due to rounding errors
        return count-1;
        }
/**
 * Sets the mean.
 */
    public void setMean(double mean) {
        if (mean != this.mean) {
            this.mean = mean;
            if (mean == -1.0) return; // not defined
            if (mean < SWITCH_MEAN) {
                this.cached_g = Math.exp(-mean);
                }
            else {
                this.cached_sq = Math.sqrt(2.0*mean);
                this.cached_alxm = Math.log(mean);
                this.cached_g = mean*cached_alxm - logGamma(mean + 1.0);
                }
            }
        }

/**
 * Returns a String representation of the receiver.
 */
    public String toString() {
        return this.getClass().getName()+"("+mean+")";
        }
    }