@Override AbstractDistribution create(Random random) { return new Gamma(0.1, 0.1, random); } },
@Test public void testNextDouble() { double[] z = new double[100000]; Random gen = RandomUtils.getRandom(); for (double alpha : new double[]{1, 2, 10, 0.1, 0.01, 100}) { Gamma g = new Gamma(alpha, 1, gen); for (int i = 0; i < z.length; i++) { z[i] = g.nextDouble(); } Arrays.sort(z); // verify that empirical CDF matches theoretical one pretty closely for (double q : seq(0.01, 1, 0.01)) { double p = z[(int) (q * z.length)]; assertEquals(q, g.cdf(p), 0.01); } } }
/** * Constructs a Gamma distribution with a given shape (alpha) and rate (beta). * * @param alpha The shape parameter. * @param rate The rate parameter. * @param randomGenerator The random number generator that generates bits for us. * @throws IllegalArgumentException if <tt>alpha <= 0.0 || alpha <= 0.0</tt>. */ public Gamma(double alpha, double rate, Random randomGenerator) { this.alpha = alpha; this.rate = rate; setRandomGenerator(randomGenerator); }
private static void checkGammaCdf(double alpha, double beta, double... values) { Gamma g = new Gamma(alpha, beta, RandomUtils.getRandom()); int i = 0; for (double x : seq(0, 2 * alpha, 2 * alpha / 10)) { assertEquals(String.format(Locale.ENGLISH, "alpha=%.2f, i=%d, x=%.2f", alpha, i, x), values[i], g.cdf(x), 1.0e-7); i++; } }
@Test public void testPdf() { Random gen = RandomUtils.getRandom(); for (double alpha : new double[]{0.01, 0.1, 1, 2, 10, 100}) { for (double beta : new double[]{0.1, 1, 2, 100}) { Gamma g1 = new Gamma(alpha, beta, gen); for (double x : seq(0, 0.99, 0.1)) { double p = Math.pow(beta, alpha) * Math.pow(x, alpha - 1) * Math.exp(-beta * x - org.apache.mahout.math.jet.stat.Gamma.logGamma(alpha)); assertEquals(String.format(Locale.ENGLISH, "alpha=%.2f, beta=%.2f, x=%.2f\n", alpha, beta, x), p, g1.pdf(x), 1.0e-9); } } } } }
@Test(timeout=50000) public void testTimesSparseEfficiency() { Random raw = RandomUtils.getRandom(); Gamma gen = new Gamma(0.1, 0.1, raw); int[] values = new int[1000]; for (int k = 0; k < 1000; k++) { int j = (int) Math.min(1000, gen.nextDouble()); values[j]++; int[] values = new int[1000]; for (int k = 0; k < 1000; k++) { int j = (int) Math.min(1000, gen.nextDouble()); values[j]++;
/** Returns a random number from the distribution. */ @Override public double nextDouble() { return nextDouble(alpha, rate); }
/** Returns the probability distribution function. * @param x Where to compute the density function. * * @return The value of the gamma density at x. */ @Override public double pdf(double x) { if (x < 0) { throw new IllegalArgumentException(); } if (x == 0) { if (alpha == 1.0) { return rate; } else if (alpha < 1) { return Double.POSITIVE_INFINITY; } else { return 0; } } if (alpha == 1.0) { return rate * Math.exp(-x * rate); } return rate * Math.exp((alpha - 1.0) * Math.log(x * rate) - x * rate - logGamma(alpha)); }
b = 1.0 + 0.36788794412 * alpha; // Step 1 while (true) { double p = b * randomDouble(); if (p <= 1.0) { // Step 2. Case gds <= 1 gds = Math.exp(Math.log(p) / alpha); if (Math.log(randomDouble()) <= -gds) { return gds / rate; if (Math.log(randomDouble()) <= (alpha - 1.0) * Math.log(gds)) { return gds / rate; double v1; do { v1 = 2.0 * randomDouble() - 1.0; double v2 = 2.0 * randomDouble() - 1.0; v12 = v1 * v1 + v2 * v2; } while (v12 > 1.0); double u = randomDouble(); if (d * u <= t * t * t) { return gds / rate; double e; do { e = -Math.log(randomDouble()); u = randomDouble(); u = u + u - 1.0; sign_u = u > 0 ? 1.0 : -1.0;
@Test(timeout=50000) public void testTimesDenseEfficiency() { Random raw = RandomUtils.getRandom(); Gamma gen = new Gamma(0.1, 0.1, raw); int[] values = new int[1000]; for (int k = 0; k < 1000; k++) { int j = (int) Math.min(1000, gen.nextDouble()); values[j]++;
Gamma g1 = new Gamma(1, beta, gen); Gamma g2 = new Gamma(1, 1, gen); for (double x : seq(0, 0.99, 0.1)) { assertEquals(String.format(Locale.ENGLISH, "Rate invariance: x = %.4f, alpha = 1, beta = %.1f", x, beta), 1 - Math.exp(-x * beta), g1.cdf(x), 1.0e-9); assertEquals(String.format(Locale.ENGLISH, "Rate invariance: x = %.4f, alpha = 1, beta = %.1f", x, beta), g2.cdf(beta * x), g1.cdf(x), 1.0e-9); Gamma g = new Gamma(alpha, 1, gen); for (double beta : new double[]{0.1, 1, 2, 100}) { Gamma g1 = new Gamma(alpha, beta, gen); for (double x : seq(0, 0.9999, 0.001)) { assertEquals( String.format(Locale.ENGLISH, "Rate invariance: x = %.4f, alpha = %.2f, beta = %.1f", x, alpha, beta), g.cdf(x * beta), g1.cdf(x), 0);
/** * Returns a sample from this distribution. The value returned will * be the number of negative samples required before achieving r * positive samples. Each successive sample is taken independently * from a Bernouli process with probability p of success. * * The algorithm used is taken from J.H. Ahrens, U. Dieter (1974): * Computer methods for sampling from gamma, beta, Poisson and * binomial distributions, Computing 12, 223--246. * * This algorithm is essentially the same as described at * http://en.wikipedia.org/wiki/Negative_binomial_distribution#Gamma.E2.80.93Poisson_mixture * except that the notion of positive and negative outcomes is uniformly * inverted. Because the inversion is complete and consistent, this * definition is effectively identical to that defined on wikipedia. */ public int nextInt(int r, double p) { return this.poisson.nextInt(gamma.nextDouble(r, p / (1.0 - p))); }
@Test public void testPdf() { Random gen = RandomUtils.getRandom(); for (double alpha : new double[]{0.01, 0.1, 1, 2, 10, 100}) { for (double beta : new double[]{0.1, 1, 2, 100}) { Gamma g1 = new Gamma(alpha, beta, gen); for (double x : seq(0, 0.99, 0.1)) { double p = Math.pow(beta, alpha) * Math.pow(x, alpha - 1) * Math.exp(-beta * x - org.apache.mahout.math.jet.stat.Gamma.logGamma(alpha)); assertEquals(String.format(Locale.ENGLISH, "alpha=%.2f, beta=%.2f, x=%.2f\n", alpha, beta, x), p, g1.pdf(x), 1.0e-9); } } } } }
/** Returns the probability distribution function. * @param x Where to compute the density function. * * @return The value of the gamma density at x. */ @Override public double pdf(double x) { if (x < 0) { throw new IllegalArgumentException(); } if (x == 0) { if (alpha == 1.0) { return rate; } else if (alpha < 1) { return Double.POSITIVE_INFINITY; } else { return 0; } } if (alpha == 1.0) { return rate * Math.exp(-x * rate); } return rate * Math.exp((alpha - 1.0) * Math.log(x * rate) - x * rate - logGamma(alpha)); }
b = 1.0 + 0.36788794412 * alpha; // Step 1 while (true) { double p = b * randomDouble(); if (p <= 1.0) { // Step 2. Case gds <= 1 gds = Math.exp(Math.log(p) / alpha); if (Math.log(randomDouble()) <= -gds) { return gds / rate; if (Math.log(randomDouble()) <= (alpha - 1.0) * Math.log(gds)) { return gds / rate; double v1; do { v1 = 2.0 * randomDouble() - 1.0; double v2 = 2.0 * randomDouble() - 1.0; v12 = v1 * v1 + v2 * v2; } while (v12 > 1.0); double u = randomDouble(); if (d * u <= t * t * t) { return gds / rate; double e; do { e = -Math.log(randomDouble()); u = randomDouble(); u = u + u - 1.0; sign_u = u > 0 ? 1.0 : -1.0;
/** * Constructs a Negative Binomial distribution which describes the probability of getting * a particular number of negative trials (k) before getting a fixed number of positive * trials (r) where each positive trial has probability (p) of being successful. * * @param r the required number of positive trials. * @param p the probability of success. * @param randomGenerator a uniform random number generator. */ public NegativeBinomial(int r, double p, Random randomGenerator) { setRandomGenerator(randomGenerator); this.r = r; this.p = p; this.gamma = new Gamma(r, 1, randomGenerator); this.poisson = new Poisson(0.0, randomGenerator); }
@Test(timeout=50000) public void testTimesOtherSparseEfficiency() { Random raw = RandomUtils.getRandom(); Gamma gen = new Gamma(0.1, 0.1, raw); // build a sequential sparse matrix and a diagonal matrix and multiply them Matrix x = new SparseRowMatrix(1000, 2000, false); for (int i = 0; i < 1000; i++) { int[] values = new int[1000]; for (int k = 0; k < 1000; k++) { int j = (int) Math.min(1000, gen.nextDouble()); values[j]++; } for (int j = 0; j < 1000; j++) { if (values[j] > 0) { x.set(i, j, values[j]); } } } Vector d = new DenseVector(2000).assign(Functions.random()); Matrix y = new DiagonalMatrix(d); long t0 = System.nanoTime(); Matrix z = x.times(y); double elapsedTime = (System.nanoTime() - t0) * 1e-6; System.out.printf("done in %.1f ms\n", elapsedTime); for (MatrixSlice row : z) { for (Vector.Element element : row.nonZeroes()) { assertEquals(x.get(row.index(), element.index()) * d.get(element.index()), element.get(), 1e-12); } } }
@Test public void testNextDouble() { double[] z = new double[100000]; Random gen = RandomUtils.getRandom(); for (double alpha : new double[]{1, 2, 10, 0.1, 0.01, 100}) { Gamma g = new Gamma(alpha, 1, gen); for (int i = 0; i < z.length; i++) { z[i] = g.nextDouble(); } Arrays.sort(z); // verify that empirical CDF matches theoretical one pretty closely for (double q : seq(0.01, 1, 0.01)) { double p = z[(int) (q * z.length)]; assertEquals(q, g.cdf(p), 0.01); } } }
private static void checkGammaCdf(double alpha, double beta, double... values) { Gamma g = new Gamma(alpha, beta, RandomUtils.getRandom()); int i = 0; for (double x : seq(0, 2 * alpha, 2 * alpha / 10)) { assertEquals(String.format(Locale.ENGLISH, "alpha=%.2f, i=%d, x=%.2f", alpha, i, x), values[i], g.cdf(x), 1.0e-7); i++; } }
/** Returns a random number from the distribution. */ @Override public double nextDouble() { return nextDouble(alpha, rate); }