Java Gamma gammaRand(double a)

Description

Gamma random generator.

License

Open Source License

Declaration

public static double gammaRand(double a) 

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 c1 = .01;
    private static final double c2 = .222222;
    private static final double c3 = .32;
    private static final double c4 = .4;
    private static final double c5 = 1.24;
    private static final long MASK = 4294967295L;
    private static long seedi = 123456789L, seedj = 362436069L;

    /**/*from  w  ww .j a v a  2  s.  c  o m*/
     * Gamma random generator.
     */
    public static double gammaRand(double a) {

        double e, x, u0, u1, u2, v, w, c, c1, c2, c3, c4, c5;
        boolean done;

        e = Math.exp(1.0);

        if (a < 1.0) {

            /* Ahrens and Dieter algorithm */
            done = false;
            c = (a + e) / e;

            do {
                u0 = uniformRand();
                u1 = uniformRand();
                v = c * u0;

                if (v <= 1.0) {
                    x = Math.exp(Math.log(v) / a);

                    if (u1 <= Math.exp(-x)) {
                        done = true;
                    }
                } else {
                    x = -Math.log((c - v) / a);

                    if ((x > 0.0)
                            && (u1 < Math.exp((a - 1.0) * Math.log(x)))) {
                        done = true;
                    }
                }
            } while (!done);
        } else if (a == 1.0) {
            x = -Math.log(uniformRand());
        } else {

            /* Cheng and Feast algorithm */
            c1 = a - 1.0;
            c2 = (a - 1.0 / (6.0 * a)) / c1;
            c3 = 2.0 / c1;
            c4 = 2.0 / (a - 1.0) + 2.0;
            c5 = 1.0 / Math.sqrt(a);

            do {
                do {
                    u1 = uniformRand();
                    u2 = uniformRand();

                    if (a > 2.5) {
                        u1 = u2 + c5 * (1.0 - 1.86 * u1);
                    }
                } while ((u1 <= 0.0) || (u1 >= 1.0));

                w = c2 * u2 / u1;
            } while ((c3 * u1 + w + 1.0 / w) > c4
                    && (c3 * Math.log(u1) - Math.log(w) + w) > 1.0);

            x = c1 * w;
        }

        return (x);
    }

    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));
    }
}

Related

1. gammaFwd(double lc)
2. gammaln(double x)
3. gammaOfArgOn2Plus1(int d)
4. gammaOfArgOn2Plus1IncludeDivisor(int d, double divisor)
5. gammaPdf(double x, double a)