Java Binomial Coefficients binomialRand(int n, double pp)

Description

Binomial random generator from Numerical Recipes

License

Open Source License

Declaration

public static int binomialRand(int n, double pp) 

Method Source Code

//package com.java2s;
// it under the terms of the GNU General Public License as published by      //

public class Main {
    private static final double[] cof = { 76.18009172947146,
            -86.50532032941677, 24.01409824083091, -1.231739572450155,
            0.1208650973866179e-2, -0.5395239384953e-5 };
    private static final long MASK = 4294967295L;
    private static long seedi = 123456789L, seedj = 362436069L;

    /**/*from  w w  w .  ja  v a2  s  .  co  m*/
     * Binomial random generator from Numerical Recipes
     */
    public static int binomialRand(int n, double pp) {

        int j, k;
        double am, em, g, p, sq, t, y;
        double pc, plog, pclog, en;

        p = (pp <= 0.5) ? pp : 1.0 - pp;
        am = n * p;

        if (p == 0.0) {
            k = 0;
        } else if (p == 1.0) {
            k = n;
        } else if (n < 50) {
            k = 0;

            for (j = 0; j < n; j++) {
                if (uniformRand() < p) {
                    k++;
                }
            }
        } else if (am < 1.0) {
            g = Math.exp(-am);
            t = 1.0;
            k = -1;

            do {
                k++;

                t *= uniformRand();
            } while (t > g);

            if (k > n) {
                k = n;
            }
        } else {
            en = n;
            g = lngamma(en + 1.0);
            pc = 1.0 - p;
            plog = Math.log(p);
            pclog = Math.log(pc);
            sq = Math.sqrt(2.0 * am * pc);

            do {
                do {
                    y = Math.tan(Math.PI * uniformRand());
                    em = sq * y + am;
                } while ((em < 0.0) || (em >= en + 1.0));

                em = Math.floor(em);
                t = 1.2
                        * sq
                        * (1.0 + y * y)
                        * Math.exp(g - lngamma(em + 1.0)
                                - lngamma(en - em + 1.0) + em * plog
                                + (en - em) * pclog);
            } while (uniformRand() > t);

            k = (int) em;
        }

        if (p != pp) {
            k = n - k;
        }

        return (k);
    }

    public static double uniformRand() {

        seedi = (seedi * 69069 + 23606797) & MASK;
        seedj ^= (seedj << 13) & MASK;
        seedj ^= (seedj >> 17) & MASK;
        seedj ^= (seedj << 5) & MASK;

        return (((seedi + seedj) & MASK) * Math.pow(2.0, -32.0));
    }

    /**
     * This is a more literal (that is, exact) copy of the log gamma method from
     * Numerical Recipes than the following one.  It was created by cutting and
     * pasting from the PDF version of the book and then converting C syntax to
     * Java. </p> The static double array above goes with this. </p> Converted
     * to Java by Frank Wimberly
     *
     * @param xx
     * @return the value ln[?(xx)] for xx > 0
     */
    public static double lngamma(double xx) {
        //Returns the value ln[?(xx)] for xx > 0.

        if (xx <= 0)
            return Double.NaN;

        //Internal arithmetic will be done in double precision, a nicety that you can omit if ?ve-?gure
        //accuracy is good enough.
        double x, y, tmp, ser;

        int j;
        y = x = xx;
        tmp = x + 5.5;
        tmp -= (x + 0.5) * Math.log(tmp);
        ser = 1.000000000190015;
        for (j = 0; j <= 5; j++) {
            ser += cof[j] / ++y;
        }
        return -tmp + Math.log(2.5066282746310005 * ser / x);
    }
}

Related

1. binomialCoefficient(long n, long k)
2. binomialCoefficientDouble(final int n, final int k)
3. binomialCoefficientLog(final int n, final int k)
4. binomialCoefficients(int n, int k)
5. binomialPmf(int k, int n, double p)