GammaDistributionTest.java

package test.dr.distibutions;

import dr.math.distributions.GammaDistribution;
import dr.math.functionEval.GammaFunction;
import junit.framework.TestCase;
import org.apache.commons.math.ConvergenceException;
import org.apache.commons.math.FunctionEvaluationException;
import org.apache.commons.math.analysis.UnivariateRealFunction;
import org.apache.commons.math.analysis.integration.RombergIntegrator;
import org.apache.commons.math.analysis.integration.UnivariateRealIntegrator;

import java.util.Random;


public class GammaDistributionTest extends TestCase{

	/**
	 * This test stochastically draws gamma
	 * variates and compares the coded pdf 
	 * with the actual pdf.  
	 * The tolerance is required to be at most 1e-10.
	 */

    static double mypdf(double value, double shape, double scale) {
        return Math.exp((shape-1) * Math.log(value) - value/scale - GammaFunction.logGamma(shape) - shape * Math.log(scale) );
    }

	public void testPdf() throws FunctionEvaluationException  {

        final int numberOfTests = 100;
        double totErr = 0;
        double ptotErr = 0; int np = 0;
        double qtotErr = 0;

        Random random = new Random(37);

        for(int i = 0; i < numberOfTests; i++){
            final double mean = .01 + (3-0.01) * random.nextDouble();
            final double var = .01 + (3-0.01) * random.nextDouble();

            final double scale = var / mean;
            final double shape = mean / scale;

            final GammaDistribution gamma = new GammaDistribution(shape,scale);

            final double value = gamma.nextGamma();

            final double mypdf = mypdf(value, shape, scale);
            final double pdf = gamma.pdf(value);
            if ( Double.isInfinite(mypdf) && Double.isInfinite(pdf)) {
                continue;
            }

            assertFalse(Double.isNaN(mypdf));
            assertFalse(Double.isNaN(pdf));

            totErr +=  mypdf != 0 ? Math.abs((pdf - mypdf)/mypdf) : pdf;

            assertFalse("nan", Double.isNaN(totErr));
            //assertEquals("" + shape + "," + scale + "," + value, mypdf,gamma.pdf(value),1e-10);

            final double cdf = gamma.cdf(value);
            UnivariateRealFunction f = new UnivariateRealFunction() {
                public double value(double v) throws FunctionEvaluationException {
                    return mypdf(v, shape, scale);
                }
            };
            final UnivariateRealIntegrator integrator = new RombergIntegrator();
            integrator.setAbsoluteAccuracy(1e-14);
            integrator.setMaximalIterationCount(16);  // fail if it takes too much time

            double x;
            try {
                x = integrator.integrate(f, 0, value);
                ptotErr += cdf != 0.0 ? Math.abs(x-cdf)/cdf : x;
                np += 1;
                //assertTrue("" + shape + "," + scale + "," + value + " " + Math.abs(x-cdf)/x + "> 1e-6", Math.abs(1-cdf/x) < 1e-6);

                //System.out.println(shape + ","  + scale + " " + value);
            } catch( ConvergenceException e ) {
                 // can't integrate , skip test
              //  System.out.println(shape + ","  + scale + " skipped");
            }

            final double q = gamma.quantile(cdf);
            qtotErr += q != 0 ? Math.abs(q-value)/q : value;
           // assertEquals("" + shape + "," + scale + "," + value + " " + Math.abs(q-value)/value, q, value, 1e-6);
        }
        //System.out.println( !Double.isNaN(totErr) );
       // System.out.println(totErr);
        // bad test, but I can't find a good threshold that works for all individual cases 
        assertTrue("failed " + totErr/numberOfTests, totErr/numberOfTests < 1e-7);
        assertTrue("failed " + ptotErr/np, np > 0 ? (ptotErr/np < 1e-5) : true);
        assertTrue("failed " + qtotErr/numberOfTests , qtotErr/numberOfTests < 1e-7);
	}
}