Java tutorial
package xyz.lejon.sampling; import java.util.ArrayList; import java.util.Collections; import java.util.List; import java.util.concurrent.ThreadLocalRandom; //This code is from the beast-mcmc package with the following copyright /* * GammaDistribution.java * * Copyright (C) 2002-2006 Alexei Drummond and Andrew Rambaut * * This file is part of BEAST. * See the NOTICE file distributed with this work for additional * information regarding copyright ownership and licensing. * * BEAST is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as * published by the Free Software Foundation; either version 2 * of the License, or (at your option) any later version. * * BEAST is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU Lesser General Public License for more details. * * You should have received a copy of the GNU Lesser General Public * License along with BEAST; if not, write to the * Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, * Boston, MA 02110-1301 USA */ import xyz.lejon.math.GammaFunction; //This class consists of code from beast-mcmc project with the following Copyright /* * GammaDistribution.java * * Copyright (c) 2002-2013 Alexei Drummond, Andrew Rambaut and Marc Suchard * * This file is part of BEAST. * See the NOTICE file distributed with this work for additional * information regarding copyright ownership and licensing. * * BEAST is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as * published by the Free Software Foundation; either version 2 * of the License, or (at your option) any later version. * * BEAST is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU Lesser General Public License for more details. * * You should have received a copy of the GNU Lesser General Public * License along with BEAST; if not, write to the * Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, * Boston, MA 02110-1301 USA */ /** * gamma distribution. * <p/> * (Parameters: shape, scale; mean: scale*shape; variance: scale^2*shape) * * @author Korbinian Strimmer * @author Gerton Lunter * @version $Id: GammaDistribution.java,v 1.9 2006/03/30 11:12:47 rambaut Exp $ */ public class GammaDistribution { // // Public stuff // /** * Constructor */ public GammaDistribution(double shape, double scale) { this.shape = shape; this.scale = scale; this.samples = 0; } public double getShape() { return shape; } public void setShape(double value) { shape = value; } public double getScale() { return scale; } public void setScale(double value) { scale = value; } public double pdf(double x) { return pdf(x, shape, scale); } public double logPdf(double x) { return logPdf(x, shape, scale); } public double cdf(double x) { return cdf(x, shape, scale); } public double quantile(double y) { return quantile(y, shape, scale); } public double mean() { return mean(shape, scale); } public double variance() { return variance(shape, scale); } public double nextGamma() { return nextGamma(shape, scale); } /** * probability density function of the Gamma distribution * * @param x argument * @param shape shape parameter * @param scale scale parameter * @return pdf value */ public static double pdf(double x, double shape, double scale) { // return Math.pow(scale,-shape)*Math.pow(x, shape-1.0)/ // Math.exp(x/scale + GammaFunction.lnGamma(shape)); if (x < 0) throw new IllegalArgumentException(); if (x == 0) { if (shape == 1.0) return 1.0 / scale; else return 0.0; } final double xs = x / scale; if (shape == 1.0) { return Math.exp(-xs) / scale; } final double a = Math.exp((shape - 1.0) * Math.log(xs) - xs - GammaFunction.lnGamma(shape)); return a / scale; } /** * the natural log of the probability density function of the distribution * * @param x argument * @param shape shape parameter * @param scale scale parameter * @return log pdf value */ public static double logPdf(double x, double shape, double scale) { // double a = Math.pow(scale,-shape) * Math.pow(x, shape-1.0); // double b = x/scale + GammaFunction.lnGamma(shape); // return Math.log(a) - b; // AR - changed this to return -ve inf instead of throwing an // exception... This makes things // much easier when using this to calculate log likelihoods. // if (x < 0) throw new IllegalArgumentException(); if (x < 0) return Double.NEGATIVE_INFINITY; if (x == 0) { if (shape == 1.0) return Math.log(1.0 / scale); else return Double.NEGATIVE_INFINITY; } if (shape == 1.0) { return (-x / scale) - Math.log(scale); } if (shape == 0.0) // uninformative return -Math.log(x); return ((shape - 1.0) * Math.log(x / scale) - x / scale - GammaFunction.lnGamma(shape)) - Math.log(scale); } /** * cumulative density function of the Gamma distribution * * @param x argument * @param shape shape parameter * @param scale scale parameter * @return cdf value */ public static double cdf(double x, double shape, double scale) { return GammaFunction.incompleteGammaP(shape, x / scale); } /** * quantile (inverse cumulative density function) of the Gamma distribution * * @param y argument * @param shape shape parameter * @param scale scale parameter * @return icdf value */ public static double quantile(double y, double shape, double scale) { return 0.5 * scale * pointChi2(y, 2.0 * shape); } /** * mean of the Gamma distribution * * @param shape shape parameter * @param scale scale parameter * @return mean */ public static double mean(double shape, double scale) { return scale * shape; } /** * variance of the Gamma distribution. * * @param shape shape parameter * @param scale scale parameter * @return variance */ public static double variance(double shape, double scale) { return scale * scale * shape; } /** * sample from the Gamma distribution. This could be calculated using * quantile, but current algorithm is faster. * * @param shape shape parameter * @param scale scale parameter * @return sample */ public static double nextGamma(double shape, double scale) { return nextGamma(shape, scale, false); } public static double nextGamma(double shape, double scale, boolean slowCode) { double sample = 0.0; if (shape < 0.00001) { if (shape < 0) { System.out.println("Negative shape parameter"); throw new IllegalArgumentException("Negative shape parameter"); } /* * special case: shape==0.0 is an improper distribution; but * sampling works if very small values are ignored (v. large ones * don't happen) This is useful e.g. for sampling from the truncated * Gamma(0,x)-distribution. */ double minimum = 1.0e-20; double maximum = 50; double normalizingConstant = Math.log(maximum) - Math.log(minimum); // Draw from 1/x (with boundaries), and shape by exp(-x) do { sample = Math .exp(Math.log(minimum) + normalizingConstant * ThreadLocalRandom.current().nextDouble()); } while (Math.exp(-sample) < ThreadLocalRandom.current().nextDouble()); // This distribution is actually scale-free, so multiplying by // 'scale' is not necessary return sample; } if (slowCode && Math.floor(shape) == shape && shape > 4.0) { for (int i = 0; i < shape; i++) sample += -Math.log(ThreadLocalRandom.current().nextDouble()); return sample * scale; } else { // Fast special cases if (shape == 1.0) { return -Math.log(ThreadLocalRandom.current().nextDouble()) * scale; } if (shape == 2.0) { return -Math .log(ThreadLocalRandom.current().nextDouble() * ThreadLocalRandom.current().nextDouble()) * scale; } if (shape == 3.0) { return -Math.log(ThreadLocalRandom.current().nextDouble() * ThreadLocalRandom.current().nextDouble() * ThreadLocalRandom.current().nextDouble()) * scale; } if (shape == 4.0) { return -Math.log(ThreadLocalRandom.current().nextDouble() * ThreadLocalRandom.current().nextDouble() * ThreadLocalRandom.current().nextDouble() * ThreadLocalRandom.current().nextDouble()) * scale; } } // general case do { try { sample = quantile(ThreadLocalRandom.current().nextDouble(), shape, scale); } catch (IllegalArgumentException e) { // random doubles do go outside the permissible range 0.000002 < // q < 0.999998 sample = 0.0; } } while (sample == 0.0); return sample; } /** * Sample from the gamma distribution, modified by a factor exp( * -(x*bias)^-1 ), i.e. from x^(shape - 1) exp(-x/scale) exp(-1/(bias*x)) * <p/> * Works by rejection sampling, using a shifted ordinary Gamma as proposal * distribution. * <p/> * (For an older, less efficient algorithm that breaks down for shape <= 1.0, * see revision 287) */ public static double nextExpGamma(double shape, double scale, double bias) { return nextExpGamma(shape, scale, bias, false); } public static double nextExpGamma(double shape, double scale, double bias, boolean slowCode) { double sample; double accept; int iters = 0; if (slowCode) { // for testing purposes -- this can get stuck completely for small bias parameters do { sample = nextGamma(shape, scale); accept = Math.exp(-1.0 / (bias * sample)); } while (ThreadLocalRandom.current().nextDouble() > accept); } else { if (shape < 0) { //return scale / (nextExpGamma(-shape, scale, bias) * bias); return 1.0 / nextExpGamma(-shape, bias, scale); } if (shape == 0) { // sample from the restriction to x >= median, and sample the other half by using // the self-similarity transformation x -> scale / (bias*x) double median = Math.sqrt(scale / bias); double rejection_mode = 1.0 / bias; // mode of rejection distribution, x^-1 e^(-1/bias x) if (rejection_mode < median) rejection_mode = median; double rejection_norm = (1.0 / rejection_mode) * Math.exp(-1.0 / (bias * rejection_mode)); do { sample = nextGamma(1.0, scale) + median; accept = (1.0 / sample) * Math.exp(-1.0 / (sample * bias)) / rejection_norm; iters += 1; } while (ThreadLocalRandom.current().nextDouble() > accept && iters < 10000); if (iters == 10000) { System.out.println( "Severe Warning: nextExpGamma (shape=0) failed to generate a sample - returning bogus value!"); } if (ThreadLocalRandom.current().nextDouble() > 0.5) { sample = scale / (bias * sample); } return sample; } // Current algorithm works for all shape parameters > 0, but it becomes very inefficient // for v. small shape parameters. if (shape <= 0.0) { System.out.println("nextExpGamma: Illegal argument (shape parameter is must be positive)"); throw new IllegalArgumentException(""); } // the function -1/(bias*x) is bounded by x/(bias x0^2) + C, with C = -2/(bias*x0), so that these functions // coincide precisely when x=x0. This gives the scale parameter of the majorating Gamma distribution // // First calculate the value that maximizes the acceptance probability at the mean of the proposal distribution double x0 = (shape * scale + Math.sqrt(4 * scale / bias + shape * shape * scale * scale)) / 2.0; // calculate the scale parameter of the majorating Gamma distribution double majorandScale = 1.0 / ((1.0 / scale) - 1.0 / (bias * x0 * x0)); // now do rejection sampling do { sample = nextGamma(shape, majorandScale); accept = Math.exp(-(sample / x0 - 1) * (sample / x0 - 1) / (bias * sample)); iters += 1; } while (ThreadLocalRandom.current().nextDouble() > accept && iters < 10000); if (accept > 1.0) { System.out.println("PROBLEM!! This should be impossible!! Contact the authors."); } if (majorandScale < 0.0) { System.out.println("PROBLEM!! This should be impossible too!! Contact the authors."); } if (iters == 10000) { System.out.println( "Severe Warning: nextExpGamma failed to generate a sample - returning bogus value!"); } } return sample; } // Private private static double pointChi2(double prob, double v) { // Returns z so that Prob{x<z}=prob where x is Chi2 distributed with df // = v // RATNEST FORTRAN by // Best DJ & Roberts DE (1975) The percentage points of the // Chi2 distribution. Applied Statistics 24: 385-388. (AS91) final double e = 0.5e-6, aa = 0.6931471805, p = prob; double ch, a, q, p1, p2, t, x, b, s1, s2, s3, s4, s5, s6; double epsi = .01; if (p < 0.000002 || p > 1 - 0.000002) { epsi = .000001; } // if (p < 0.000002 || p > 0.999998 || v <= 0) { // throw new IllegalArgumentException("Arguments out of range p" + p + " v " + v); // } double g = GammaFunction.lnGamma(v / 2); double xx = v / 2; double c = xx - 1; if (v < -1.24 * Math.log(p)) { ch = Math.pow((p * xx * Math.exp(g + xx * aa)), 1 / xx); if (ch - e < 0) { return ch; } } else { if (v > 0.32) { x = jdistlib.Normal.quantile(p, 0, 1, true, false); //x = NormalDistribution.quantile(p, 0, 1); p1 = 0.222222 / v; ch = v * Math.pow((x * Math.sqrt(p1) + 1 - p1), 3.0); if (ch > 2.2 * v + 6) { ch = -2 * (Math.log(1 - p) - c * Math.log(.5 * ch) + g); } } else { ch = 0.4; a = Math.log(1 - p); do { q = ch; p1 = 1 + ch * (4.67 + ch); p2 = ch * (6.73 + ch * (6.66 + ch)); t = -0.5 + (4.67 + 2 * ch) / p1 - (6.73 + ch * (13.32 + 3 * ch)) / p2; ch -= (1 - Math.exp(a + g + .5 * ch + c * aa) * p2 / p1) / t; } while (Math.abs(q / ch - 1) - epsi > 0); } } do { q = ch; p1 = 0.5 * ch; if ((t = GammaFunction.incompleteGammaP(xx, p1, g)) < 0) { throw new IllegalArgumentException("Arguments out of range: t < 0"); } p2 = p - t; t = p2 * Math.exp(xx * aa + g + p1 - c * Math.log(ch)); b = t / ch; a = 0.5 * t - b * c; s1 = (210 + a * (140 + a * (105 + a * (84 + a * (70 + 60 * a))))) / 420; s2 = (420 + a * (735 + a * (966 + a * (1141 + 1278 * a)))) / 2520; s3 = (210 + a * (462 + a * (707 + 932 * a))) / 2520; s4 = (252 + a * (672 + 1182 * a) + c * (294 + a * (889 + 1740 * a))) / 5040; s5 = (84 + 264 * a + c * (175 + 606 * a)) / 2520; s6 = (120 + c * (346 + 127 * c)) / 5040; ch += t * (1 + 0.5 * t * s1 - b * c * (s1 - b * (s2 - b * (s3 - b * (s4 - b * (s5 - b * s6)))))); } while (Math.abs(q / ch - 1) > e); return (ch); } public static void main(String[] args) { testQuantile(1e-10, 0.878328435043444, 0.0013696236839573005); testQuantile(0.5, 0.878328435043444, 0.0013696236839573005); testQuantile(1.0 - 1e-10, 0.878328435043444, 0.0013696236839573005); //testQuantileCM(1e-10, 0.878328435043444, 0.0013696236839573005); //testQuantileCM(0.5, 0.878328435043444, 0.0013696236839573005); //testQuantileCM(1.0 - 1e-10, 0.878328435043444, 0.0013696236839573005); /* for (double i = 0.0125; i < 1.0; i += 0.025) { System.out.print(i + ": "); try { System.out.println(new GammaDistributionImpl(0.878328435043444, 0.0013696236839573005).inverseCumulativeProbability(i)); } catch (MathException e) { System.out.println(e.getMessage()); } }*/ // System.out // .println("K-S critical values: 1.22(10%), 1.36(5%), 1.63(1%)\n"); // // int iters = 30000; // testExpGamma(1.0, 0.01, 7, iters); // testExpGamma(1.0, 0.01, 5, iters); // testExpGamma(2.0, 0.01, 10000, iters); // testExpGamma(1.0, 0.01, 10000, iters); // testExpGamma(0.1, 0.01, 10000, iters); // testExpGamma(0.01, 0.01, 10000, iters); // // testExpGamma(2.0, 0.01, 10, iters); // testExpGamma(1.5, 0.01, 10, iters); // testExpGamma(1.0, 0.01, 10, iters); // testExpGamma(0.9, 0.01, 10, iters); // testExpGamma(0.5, 0.01, 10, iters); // testExpGamma(0.4, 0.01, 10, iters); // testExpGamma(0.3, 0.01, 10, iters); // testExpGamma(0.2, 0.01, 10, iters); // testExpGamma(0.1, 0.01, 10, iters); // // // test distributions with severe bias, where rejection sampling doesn't // // work anymore // testExpGamma2(2.0, 0.01, 1, iters, 0.112946); // testExpGamma2(2.0, 0.01, 0.1, iters, 0.328874); // testExpGamma2(2.0, 0.01, 0.01, iters, 1.01255); // testExpGamma2(1.0, 0.01, 0.0003, iters, 5.781); // testExpGamma2(4.0, 0.01, 0.0003, iters, 5.79604); // testExpGamma2(20.0, 0.01, 0.0003, iters, 5.87687); // testExpGamma2(10.0, 0.01, 0.01, iters, 1.05374); // testExpGamma2(1.0, 0.01, 0.05, iters, 0.454734); // // test the basic Gamma distribution // test(1.0, 1.0, iters); // test(2.0, 1.0, iters); // test(3.0, 1.0, iters); // test(4.0, 1.0, iters); // test(100.0, 1.0, iters); // testAddition(0.5, 1.0, 2, iters); // testAddition(0.25, 1.0, 4, iters); // testAddition(0.1, 1.0, 10, iters); // testAddition(10, 1.0, 10, iters); // testAddition(20, 1.0, 10, iters); // test(0.001, 1.0, iters); // test(1.0, 2.0, iters); // test(10.0, 1.0, iters); // test(16.0, 1.0, iters); // test(16.0, 0.1, iters); // test(100.0, 1.0, iters); // test(0.5, 1.0, iters); // test(0.5, 0.1, iters); // test(0.1, 1.0, iters); // test(0.9, 1.0, iters); // // test distributions with milder biases, and compare with results from // // simple rejection sampling // testExpGamma(2.0, 0.000001, 1000000, iters); // testExpGamma(2.0, 0.000001, 100000, iters); // testExpGamma(2.0, 0.000001, 70000, iters); // testExpGamma(10.0, 0.01, 7, iters); // testExpGamma(10.0, 0.01, 5, iters); // testExpGamma(1.0, 0.01, 100, iters); // testExpGamma(1.0, 0.01, 10, iters); // testExpGamma(1.0, 0.01, 7, iters / 3); // testExpGamma(1.0, 0.01, 5, iters / 3); // testExpGamma(1.0, 0.00001, 1000000, iters); // testExpGamma(1.0, 0.00001, 100000, iters); // testExpGamma(1.0, 0.00001, 10000, iters); // testExpGamma(1.0, 0.00001, 5000, iters / 3); /* // * this one takes some // * time // */ // testExpGamma(2.0, 1.0, 0.5, iters); // testExpGamma(2.0, 1.0, 1.0, iters); // testExpGamma(2.0, 1.0, 2.0, iters); // testExpGamma(3.0, 3.0, 2.0, iters); // testExpGamma(10.0, 3.0, 5.0, iters); // testExpGamma(1.0, 3.0, 5.0, iters); // testExpGamma(1.0, 10.0, 5.0, iters); // testExpGamma(2.0, 10.0, 5.0, iters); // // test the basic Gamma distribution // test(1.0, 1.0, iters); // test(2.0, 1.0, iters); // test(3.0, 1.0, iters); // test(4.0, 1.0, iters); // test(100.0, 1.0, iters); // testAddition(0.5, 1.0, 2, iters); // testAddition(0.25, 1.0, 4, iters); // testAddition(0.1, 1.0, 10, iters); // testAddition(10, 1.0, 10, iters); // testAddition(20, 1.0, 10, iters); // test(0.001, 1.0, iters); // test(1.0, 2.0, iters); // test(10.0, 1.0, iters); // test(16.0, 1.0, iters); // test(16.0, 0.1, iters); // test(100.0, 1.0, iters); // test(0.5, 1.0, iters); // test(0.5, 0.1, iters); // test(0.1, 1.0, iters); // test(0.9, 1.0, iters); } private static void testQuantile(double y, double shape, double scale) { long time = System.currentTimeMillis(); double value = 0; for (int i = 0; i < 1000; i++) { value = quantile(y, shape, scale); } value = quantile(y, shape, scale); long elapsed = System.currentTimeMillis() - time; System.out.println("Quantile, " + y + ", for shape=" + shape + ", scale=" + scale + " : " + value + ", time=" + elapsed + "ms"); } /* private static void testQuantileCM(double y, double shape, double scale) { long time = System.currentTimeMillis(); double value = 0; try { for (int i = 0; i < 1000; i++) { value = (new org.apache.commons.math.distribution.GammaDistributionImpl(shape, scale)).inverseCumulativeProbability(y); } value = (new org.apache.commons.math.distribution.GammaDistributionImpl(shape, scale)).inverseCumulativeProbability(y); } catch (MathException e) { e.printStackTrace(); } long elapsed = System.currentTimeMillis() - time; System.out.println("commons.maths inverseCDF, "+ y +", for shape=" + shape + ", scale=" + scale + " : " + value + ", time=" + elapsed + "ms"); }*/ /* assumes the two arrays are the same size */ private static double KolmogorovSmirnov(List<Double> l1, List<Double> l2) { int idx2 = 0; int max = 0; for (int i = 0; i < l1.size(); i++) { while (idx2 < l2.size() && l2.get(idx2) < l1.get(i)) { idx2++; } max = Math.max(max, idx2 - i); } return max / Math.sqrt(2.0 * l1.size()); } @SuppressWarnings("unused") private static void testExpGamma2(double shape, double scale, double bias, int iterations, double mean) { double s0 = 0; double s1 = 0; double s2 = 0; List<Double> fast = new ArrayList<Double>(0); for (int i = 0; i < iterations; i++) { double sample = nextExpGamma(shape, scale, bias, false); s0 += 1; s1 += sample; s2 += sample * sample; fast.add(sample); } Collections.sort(fast); double expmean = s1 / s0; double expvar = (s2 - s1 * s1 / s0) / s0; double z = (mean - expmean) / Math.sqrt(expvar / iterations); System.out.println("Equal-mean test: (shape=" + shape + " scale=" + scale + " bias=" + bias + " mean=" + expmean + " expected=" + mean + " var=" + expvar + " median=" + fast.get(iterations / 2) + "): z=" + z); } @SuppressWarnings("unused") private static void testExpGamma(double shape, double scale, double bias, int iterations) { List<Double> slow = new ArrayList<Double>(0); List<Double> fast = new ArrayList<Double>(0); long time = System.currentTimeMillis(); for (int i = 0; i < iterations; i++) { slow.add(nextExpGamma(shape, scale, bias, true)); } long slowtime = System.currentTimeMillis(); for (int i = 0; i < iterations; i++) { fast.add(nextExpGamma(shape, scale, bias, false)); } long fasttime = System.currentTimeMillis() - slowtime; slowtime -= time; Collections.sort(slow); Collections.sort(fast); System.out.println("KS test for shape=" + shape + ", bias=" + bias + " : " + KolmogorovSmirnov(slow, fast) + " and " + KolmogorovSmirnov(fast, slow) + " slow=" + slowtime + "ms, fast=" + fasttime + "ms"); } @SuppressWarnings("unused") private static void test(double shape, double scale, int iterations) { List<Double> slow = new ArrayList<Double>(0); List<Double> fast = new ArrayList<Double>(0); for (int i = 0; i < iterations; i++) { slow.add(nextGamma(shape, scale, true)); fast.add(nextGamma(shape, scale, false)); } Collections.sort(slow); Collections.sort(fast); System.out.println("KS test for shape=" + shape + " : " + KolmogorovSmirnov(slow, fast) + " and " + KolmogorovSmirnov(fast, slow)); } @SuppressWarnings("unused") private static void testAddition(double shape, double scale, int N, int iterations) { List<Double> slow = new ArrayList<Double>(0); List<Double> fast = new ArrayList<Double>(0); List<Double> test = new ArrayList<Double>(0); for (int i = 0; i < iterations; i++) { double s = 0.0; for (int j = 0; j < N; j++) { s += nextGamma(shape, scale, true); } slow.add(s); s = 0.0; for (int j = 0; j < N; j++) { s += nextGamma(shape, scale, false); } fast.add(s); test.add(nextGamma(shape * N, scale, true)); } Collections.sort(slow); Collections.sort(fast); Collections.sort(test); System.out.println("KS test for shape=" + shape + " : slow=" + KolmogorovSmirnov(slow, test) + " & " + KolmogorovSmirnov(test, slow) + "; fast=" + KolmogorovSmirnov(fast, test) + " & " + KolmogorovSmirnov(test, fast)); } protected double shape, scale; protected int samples; }