Java tutorial
/* * #%L * gitools-core * %% * Copyright (C) 2013 Universitat Pompeu Fabra - Biomedical Genomics group * %% * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3 of the * License, or (at your option) any later version. * * This program 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 General Public License for more details. * * You should have received a copy of the GNU General Public * License along with this program. If not, see * <http://www.gnu.org/licenses/gpl-3.0.html>. * #L% */ package org.gitools.analysis.stats.test; import org.apache.commons.math3.distribution.BinomialDistribution; import org.apache.commons.math3.distribution.NormalDistribution; import org.apache.commons.math3.distribution.PoissonDistribution; import org.gitools.analysis.stats.test.results.BinomialResult; import org.gitools.analysis.stats.test.results.CommonResult; import static com.google.common.base.Predicates.notNull; import static com.google.common.collect.Iterables.filter; public class BinomialTest extends AbstractEnrichmentTest { private static final int exactSizeLimit = 100000; public enum AproximationMode { onlyExact, onlyNormal, onlyPoisson, automatic } private abstract class BinomialAproximation { public abstract CommonResult getResult(int observed, int n, double p, double expectedMean, double expectedStdev, double expectedVar); } private final AproximationMode aproxMode; private double p; private BinomialAproximation aprox; public BinomialTest(AproximationMode aproxMode) { super("binomial", BinomialResult.class); this.aproxMode = aproxMode; switch (aproxMode) { case onlyExact: this.aprox = new BinomialAproximation() { @Override public CommonResult getResult(int observed, int n, double p, double expectedMean, double expectedStdev, double expectedVar) { return resultWithExact(observed, n, p, expectedMean, expectedStdev); } }; break; case onlyNormal: this.aprox = new BinomialAproximation() { @Override public CommonResult getResult(int observed, int n, double p, double expectedMean, double expectedStdev, double expectedVar) { return resultWithNormal(observed, n, p, expectedMean, expectedStdev); } }; break; case onlyPoisson: this.aprox = new BinomialAproximation() { @Override public CommonResult getResult(int observed, int n, double p, double expectedMean, double expectedStdev, double expectedVar) { return resultWithPoisson(observed, n, p, expectedMean, expectedStdev); } }; break; case automatic: this.aprox = new BinomialAproximation() { @Override public CommonResult getResult(int observed, int n, double p, double expectedMean, double expectedStdev, double expectedVar) { if (n <= exactSizeLimit) { return resultWithExact(observed, n, p, expectedMean, expectedStdev); } else if (n >= 150 && expectedVar >= 0.9 * expectedMean) { return resultWithPoisson(observed, n, p, expectedMean, expectedStdev); } else if ((n * p >= 5) && (n * (1 - p) >= 5)) { return resultWithNormal(observed, n, p, expectedMean, expectedStdev); } else { return resultWithExact(observed, n, p, expectedMean, expectedStdev); } } }; break; } } @Override public String getName() { StringBuilder sb = new StringBuilder(); sb.append("binomial"); switch (aproxMode) { case automatic: break; case onlyExact: sb.append("-exact"); break; case onlyNormal: sb.append("-normal"); break; case onlyPoisson: sb.append("-poisson"); break; } return sb.toString(); } @Override public void processPopulation(Iterable<Double> population) { double size = 0; double observed = 0; for (Double value : population) { if (value == null) { continue; } if (value == 1.0) { observed++; } size++; } p = observed / size; } @Override public CommonResult processTest(Iterable<Double> values) { int observed = 0; int n = 0; for (Double value : filter(values, notNull())) { if (value == 1.0) { observed++; } n++; } double expectedMean = n * p; double expectedVar = n * p * (1.0 - p); double expectedStdev = Math.sqrt(expectedVar); return aprox.getResult(observed, n, p, expectedMean, expectedStdev, expectedVar); } private static CommonResult resultWithExact(int observed, int n, double p, double expectedMean, double expectedStdev) { double leftPvalue; double rightPvalue; double twoTailPvalue; //FIXME: May be it's better to return null ??? if (n == 0) { leftPvalue = rightPvalue = twoTailPvalue = 1.0; } else { BinomialDistribution distribution = new BinomialDistribution(n, p); leftPvalue = distribution.cumulativeProbability(observed); rightPvalue = observed > 0 ? distribution.cumulativeProbability(observed - 1, n) : 1.0; twoTailPvalue = leftPvalue + rightPvalue; twoTailPvalue = twoTailPvalue > 1.0 ? 1.0 : twoTailPvalue; } return new BinomialResult(BinomialResult.Distribution.BINOMIAL, n, leftPvalue, rightPvalue, twoTailPvalue, observed, expectedMean, expectedStdev, p); } private static NormalDistribution NORMAL = new NormalDistribution(); private static CommonResult resultWithNormal(int observed, int n, double p, double expectedMean, double expectedStdev) { double zscore; double leftPvalue; double rightPvalue; double twoTailPvalue; // Calculate zscore and pvalues zscore = (observed - expectedMean) / expectedStdev; leftPvalue = NORMAL.cumulativeProbability(zscore); rightPvalue = 1.0 - leftPvalue; twoTailPvalue = (zscore <= 0 ? leftPvalue : rightPvalue) * 2; twoTailPvalue = twoTailPvalue > 1.0 ? 1.0 : twoTailPvalue; return new BinomialResult(BinomialResult.Distribution.NORMAL, n, leftPvalue, rightPvalue, twoTailPvalue, observed, expectedMean, expectedStdev, p); } private static CommonResult resultWithPoisson(int observed, int n, double p, double expectedMean, double expectedStdev) { double leftPvalue; double rightPvalue; double twoTailPvalue; try { PoissonDistribution poisson = new PoissonDistribution(expectedMean); leftPvalue = poisson.cumulativeProbability(observed); rightPvalue = observed > 0 ? poisson.cumulativeProbability(observed - 1, n) : 1.0; twoTailPvalue = (observed <= expectedMean ? leftPvalue : rightPvalue) * 2; //FIXME: Review twoTailPvalue = twoTailPvalue > 1.0 ? 1.0 : twoTailPvalue; } catch (ArithmeticException e) { leftPvalue = rightPvalue = twoTailPvalue = Double.NaN; } return new BinomialResult(BinomialResult.Distribution.POISSON, n, leftPvalue, rightPvalue, twoTailPvalue, observed, expectedMean, expectedStdev, p); } }