commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From pste...@apache.org
Subject svn commit: r1405620 - in /commons/proper/math/trunk: pom.xml src/changes/changes.xml src/main/java/org/apache/commons/math3/stat/inference/GTest.java src/test/java/org/apache/commons/math3/stat/inference/GTestTest.java
Date Sun, 04 Nov 2012 19:32:36 GMT
Author: psteitz
Date: Sun Nov  4 19:32:35 2012
New Revision: 1405620

URL: http://svn.apache.org/viewvc?rev=1405620&view=rev
Log:
Added G-test statistics. JIRA: MATH-878.  Thanks to Radoslav Tsvetkov and Ted Dunning.

Added:
    commons/proper/math/trunk/src/main/java/org/apache/commons/math3/stat/inference/GTest.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math3/stat/inference/GTestTest.java
Modified:
    commons/proper/math/trunk/pom.xml
    commons/proper/math/trunk/src/changes/changes.xml

Modified: commons/proper/math/trunk/pom.xml
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/pom.xml?rev=1405620&r1=1405619&r2=1405620&view=diff
==============================================================================
--- commons/proper/math/trunk/pom.xml (original)
+++ commons/proper/math/trunk/pom.xml Sun Nov  4 19:32:35 2012
@@ -265,6 +265,9 @@
       <name>Mauro Talevi</name>
     </contributor>
     <contributor>
+      <name>Radoslav Tsvetkov</name>
+    </contributor>
+    <contributor>
       <name>Kim van der Linde</name>
     </contributor>
     <contributor>

Modified: commons/proper/math/trunk/src/changes/changes.xml
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/changes/changes.xml?rev=1405620&r1=1405619&r2=1405620&view=diff
==============================================================================
--- commons/proper/math/trunk/src/changes/changes.xml (original)
+++ commons/proper/math/trunk/src/changes/changes.xml Sun Nov  4 19:32:35 2012
@@ -52,6 +52,9 @@ If the output is not quite correct, chec
   <body>
     <release version="3.1" date="TBD" description="
 ">
+      <action dev="psteitz" type="add" issue="MATH-878" due-to="Radoslav Tsvetkov">
+        Added G-test statistics.
+      </action>
       <action dev="erans" type="add" issue="MATH-883">
         New "getSquareRoot" method in class "EigenDecomposition" (package
         "o.a.c.m.linear").

Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/stat/inference/GTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/stat/inference/GTest.java?rev=1405620&view=auto
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/stat/inference/GTest.java
(added)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/stat/inference/GTest.java
Sun Nov  4 19:32:35 2012
@@ -0,0 +1,537 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.math3.stat.inference;
+
+import org.apache.commons.math3.distribution.ChiSquaredDistribution;
+import org.apache.commons.math3.exception.DimensionMismatchException;
+import org.apache.commons.math3.exception.MaxCountExceededException;
+import org.apache.commons.math3.exception.NotPositiveException;
+import org.apache.commons.math3.exception.NotStrictlyPositiveException;
+import org.apache.commons.math3.exception.OutOfRangeException;
+import org.apache.commons.math3.exception.ZeroException;
+import org.apache.commons.math3.exception.util.LocalizedFormats;
+import org.apache.commons.math3.util.FastMath;
+import org.apache.commons.math3.util.MathArrays;
+
+/**
+ * Implements <a href="http://en.wikipedia.org/wiki/G-test">G Test</a>
+ * statistics.
+ *
+ * <p>This is known in statistical genetics as the McDonald-Kreitman test.
+ * The implementation handles both known and unknown distributions.</p>
+ *
+ * <p>Two samples tests can be used when the distribution is unknown <i>a priori</i>
+ * but provided by one sample, or when the hypothesis under test is that the two
+ * samples come from the same underlying distribution.</p>
+ *
+ * @version $Id$
+ * @since 3.1
+ */
+public class GTest {
+
+    /**
+     * Computes the <a href="http://en.wikipedia.org/wiki/G-test">G statistic
+     * for Goodness of Fit</a> comparing {@code observed} and {@code expected}
+     * frequency counts.
+     *
+     * <p>This statistic can be used to perform a G test (Log-Likelihood Ratio
+     * Test) evaluating the null hypothesis that the observed counts follow the
+     * expected distribution.</p>
+     *
+     * <p><strong>Preconditions</strong>: <ul>
+     * <li>Expected counts must all be positive. </li>
+     * <li>Observed counts must all be &ge; 0. </li>
+     * <li>The observed and expected arrays must have the same length and their
+     * common length must be at least 2. </li></ul></p>
+     *
+     * <p>If any of the preconditions are not met, a
+     * {@code MathIllegalArgumentException} is thrown.</p>
+     *
+     * <p><strong>Note:</strong>This implementation rescales the
+     * {@code expected} array if necessary to ensure that the sum of the
+     * expected and observed counts are equal.</p>
+     *
+     * @param observed array of observed frequency counts
+     * @param expected array of expected frequency counts
+     * @return G-Test statistic
+     * @throws NotPositiveException if {@code observed} has negative entries
+     * @throws NotStrictlyPositiveException if {@code expected} has entries that
+     * are not strictly positive
+     * @throws DimensionMismatchException if the array lengths do not match or
+     * are less than 2.
+     */
+    public double gValueGoodnessOfFit(final double[] expected, final long[] observed)
+            throws NotPositiveException, NotStrictlyPositiveException,
+            DimensionMismatchException {
+
+        if (expected.length < 2) {
+            throw new DimensionMismatchException(expected.length, 2);
+        }
+        if (expected.length != observed.length) {
+            throw new DimensionMismatchException(expected.length, observed.length);
+        }
+        MathArrays.checkPositive(expected);
+        MathArrays.checkNonNegative(observed);
+
+        double sumExpected = 0d;
+        double sumObserved = 0d;
+        for (int i = 0; i < observed.length; i++) {
+            sumExpected += expected[i];
+            sumObserved += observed[i];
+        }
+        double ratio = 1d;
+        boolean rescale = false;
+        if (Math.abs(sumExpected - sumObserved) > 10E-6) {
+            ratio = sumObserved / sumExpected;
+            rescale = true;
+        }
+        double sum = 0d;
+        for (int i = 0; i < observed.length; i++) {
+            final double dev = rescale ?
+                    FastMath.log((double) observed[i] / (ratio * expected[i])) :
+                        FastMath.log((double) observed[i] / expected[i]);
+            sum += ((double) observed[i]) * dev;
+        }
+        return 2d * sum;
+    }
+
+    /**
+     * Returns the <i>observed significance level</i>, or <a href=
+     * "http://www.cas.lancs.ac.uk/glossary_v1.1/hyptest.html#pvalue"> p-value</a>,
+     * associated with a G-Test for goodness of fit</a> comparing the
+     * {@code observed} frequency counts to those in the {@code expected} array.
+     *
+     * <p>The number returned is the smallest significance level at which one
+     * can reject the null hypothesis that the observed counts conform to the
+     * frequency distribution described by the expected counts.</p>
+     *
+     * <p>The probability returned is the tail probability beyond
+     * {@link #gValueGoodnessOfFit(double[], long[]) gValueGoodnessOfFit(expected, observed)}
+     * in the ChiSquare distribution with degrees of freedom one less than the
+     * common length of {@code expected} and {@code observed}.</p>
+     *
+     * <p> <strong>Preconditions</strong>: <ul>
+     * <li>Expected counts must all be positive. </li>
+     * <li>Observed counts must all be &ge; 0. </li>
+     * <li>The observed and expected arrays must have the
+     * same length and their common length must be at least 2.</li>
+     * </ul></p>
+     *
+     * <p>If any of the preconditions are not met, a
+     * {@code MathIllegalArgumentException} is thrown.</p>
+     *
+     * <p><strong>Note:</strong>This implementation rescales the
+     * {@code expected} array if necessary to ensure that the sum of the
+     *  expected and observed counts are equal.</p>
+     *
+     * @param observed array of observed frequency counts
+     * @param expected array of expected frequency counts
+     * @return p-value
+     * @throws NotPositiveException if {@code observed} has negative entries
+     * @throws NotStrictlyPositiveException if {@code expected} has entries that
+     * are not strictly positive
+     * @throws DimensionMismatchException if the array lengths do not match or
+     * are less than 2.
+     * @throws MaxCountExceededException if an error occurs computing the
+     * p-value.
+     */
+    public double gTestGoodnessOfFitPValue(final double[] expected, final long[] observed)
+            throws NotPositiveException, NotStrictlyPositiveException,
+            DimensionMismatchException, MaxCountExceededException {
+
+        final ChiSquaredDistribution distribution =
+                new ChiSquaredDistribution(expected.length - 1.0);
+        return 1.0 - distribution.cumulativeProbability(
+                gValueGoodnessOfFit(expected, observed));
+    }
+
+    /**
+     * Returns the intrinsic (Hardy-Weinberg proportions) p-Value, as described
+     * in p64-69 of McDonald, J.H. 2009. Handbook of Biological Statistics
+     * (2nd ed.). Sparky House Publishing, Baltimore, Maryland.
+     *
+     * <p> The probability returned is the tail probability beyond
+     * {@link #gValueGoodnessOfFit(double[], long[]) gValueGoodnessOfFit(expected, observed)}
+     * in the ChiSquare distribution with degrees of freedom two less than the
+     * common length of {@code expected} and {@code observed}.</p>
+     *
+     * @param observed array of observed frequency counts
+     * @param expected array of expected frequency counts
+     * @return p-value
+     * @throws NotPositiveException if {@code observed} has negative entries
+     * @throws NotStrictlyPositiveException {@code expected} has entries that are
+     * not strictly positive
+     * @throws DimensionMismatchException if the array lengths do not match or
+     * are less than 2.
+     * @throws MaxCountExceededException if an error occurs computing the
+     * p-value.
+     */
+    public double gTestGoodnessOfFitIntrinsicPValue(final double[] expected, final long[]
observed)
+            throws NotPositiveException, NotStrictlyPositiveException,
+            DimensionMismatchException, MaxCountExceededException {
+
+        final ChiSquaredDistribution distribution =
+                new ChiSquaredDistribution(expected.length - 2.0);
+        return 1.0 - distribution.cumulativeProbability(
+                gValueGoodnessOfFit(expected, observed));
+    }
+
+    /**
+     * Performs a G-Test (Log-Likelihood Ratio Test) for goodness of fit
+     * evaluating the null hypothesis that the observed counts conform to the
+     * frequency distribution described by the expected counts, with
+     * significance level {@code alpha}. Returns true iff the null
+     * hypothesis can be rejected with {@code 100 * (1 - alpha)} percent confidence.
+     *
+     * <p><strong>Example:</strong><br> To test the hypothesis that
+     * {@code observed} follows {@code expected} at the 99% level,
+     * use </p><p>
+     * {@code gTest(expected, observed, 0.01)}</p>
+     *
+     * <p>Returns true iff {@link #gTestGoodnessOfFitPValue(double[], long[])
+     *  gTestGoodnessOfFitPValue(expected, observed)} < alpha</p>
+     *
+     * <p><strong>Preconditions</strong>: <ul>
+     * <li>Expected counts must all be positive. </li>
+     * <li>Observed counts must all be &ge; 0. </li>
+     * <li>The observed and expected arrays must have the same length and their
+     * common length must be at least 2.
+     * <li> {@code 0 < alpha < 0.5} </li></ul></p>
+     *
+     * <p>If any of the preconditions are not met, a
+     * {@code MathIllegalArgumentException} is thrown.</p>
+     *
+     * <p><strong>Note:</strong>This implementation rescales the
+     * {@code expected} array if necessary to ensure that the sum of the
+     * expected and observed counts are equal.</p>
+     *
+     * @param observed array of observed frequency counts
+     * @param expected array of expected frequency counts
+     * @param alpha significance level of the test
+     * @return true iff null hypothesis can be rejected with confidence 1 -
+     * alpha
+     * @throws NotPositiveException if {@code observed} has negative entries
+     * @throws NotStrictlyPositiveException if {@code expected} has entries that
+     * are not strictly positive
+     * @throws DimensionMismatchException if the array lengths do not match or
+     * are less than 2.
+     * @throws MaxCountExceededException if an error occurs computing the
+     * p-value.
+     * @throws OutOfRangeException if alpha is not strictly greater than zero
+     * and less than or equal to 0.5
+     */
+    public boolean gTestGoodnessOfFit(final double[] expected, final long[] observed,
+            final double alpha)
+            throws NotPositiveException, NotStrictlyPositiveException,
+            DimensionMismatchException, OutOfRangeException, MaxCountExceededException {
+
+        if ((alpha <= 0) || (alpha > 0.5)) {
+            throw new OutOfRangeException(LocalizedFormats.OUT_OF_BOUND_SIGNIFICANCE_LEVEL,
+                    alpha, 0, 0.5);
+        }
+        return gTestGoodnessOfFitPValue(expected, observed) < alpha;
+    }
+
+    /**
+     * Calculates the <a href=
+     * "http://en.wikipedia.org/wiki/Entropy_%28information_theory%29">Shannon
+     * entropy</a> for 2 Dimensional Matrix.  The value returned is the entropy
+     * of the vector formed by concatenating the rows (or columns) of {@code k}
+     * to form a vector. See {@link #entropy(long[])}.
+     *
+     * @param k 2 Dimensional Matrix of long values (for ex. the counts of a
+     * trials)
+     * @return Shannon Entropy of the given Matrix
+     *
+     */
+    private double entropy(final long[][] k) {
+        double h = 0d;
+        double sum_k = 0d;
+        for (int i = 0; i < k.length; i++) {
+            for (int j = 0; j < k[i].length; j++) {
+                sum_k += (double) k[i][j];
+            }
+        }
+        for (int i = 0; i < k.length; i++) {
+            for (int j = 0; j < k[i].length; j++) {
+                if (k[i][j] != 0) {
+                    final double p_ij = (double) k[i][j] / sum_k;
+                    h += p_ij * Math.log(p_ij);
+                }
+            }
+        }
+        return -h;
+    }
+
+    /**
+     * Calculates the <a href="http://en.wikipedia.org/wiki/Entropy_%28information_theory%29">
+     * Shannon entropy</a> for a vector.  The values of {@code k} are taken to be
+     * incidence counts of the values of a random variable. What is returned is <br/>
+     * &sum;p<sub>i</sub>log(p<sub>i</sub><br/>
+     * where p<sub>i</sub> = k[i] / (sum of elements in k)
+     *
+     * @param k Vector (for ex. Row Sums of a trials)
+     * @return Shannon Entropy of the given Vector
+     *
+     */
+    private double entropy(final long[] k) {
+        double h = 0d;
+        double sum_k = 0d;
+        for (int i = 0; i < k.length; i++) {
+            sum_k += (double) k[i];
+        }
+        for (int i = 0; i < k.length; i++) {
+            if (k[i] != 0) {
+                final double p_i = (double) k[i] / sum_k;
+                h += p_i * Math.log(p_i);
+            }
+        }
+        return -h;
+    }
+
+    /**
+     * <p>Computes a G (Log-Likelihood Ratio) two sample test statistic for
+     * independence comparing frequency counts in
+     * {@code observed1} and {@codeobserved2}. The sums of frequency
+     * counts in the two samples are not required to be the same. The formula
+     * used to compute the test statistic is </p>
+     *
+     * <p>{@code 2 * totalSum * [H(rowSums) + H(colSums) - H(k)]}</p>
+     *
+     * <p> where {@code H} is the
+     * <a href="http://en.wikipedia.org/wiki/Entropy_%28information_theory%29">
+     * Shannon Entropy</a> of the random variable formed by viewing the elements
+     * of the argument array as incidence counts; <br/>
+     * {@code k} is a matrix with rows {@code [observed1, observed2]}; <br/>
+     * {@code rowSums, colSums} are the row/col sums of {@code k}; <br>
+     * and {@code totalSum} is the overall sum of all entries in {@code k}.</p>
+     *
+     * <p>This statistic can be used to perform a G test evaluating the null
+     * hypothesis that both observed counts are independent </p>
+     *
+     * <p> <strong>Preconditions</strong>: <ul>
+     * <li>Observed counts must be non-negative. </li>
+     * <li>Observed counts for a specific bin must not both be zero. </li>
+     * <li>Observed counts for a specific sample must not all be  0. </li>
+     * <li>The arrays {@code observed1} and {@code observed2} must have
+     * the same length and their common length must be at least 2. </li></ul></p>
+     *
+     * <p>If any of the preconditions are not met, a
+     * {@code MathIllegalArgumentException} is thrown.</p>
+     *
+     * @param observed1 array of observed frequency counts of the first data set
+     * @param observed2 array of observed frequency counts of the second data
+     * set
+     * @return G-Test statistic
+     * @throws DimensionMismatchException the the lengths of the arrays do not
+     * match or their common length is less than 2
+     * @throws NotPositiveException if any entry in {@code observed1} or
+     * {@code observed2} is negative
+     * @throws ZeroException if either all counts of
+     * {@code observed1} or {@code observed2} are zero, or if the count
+     * at the same index is zero for both arrays.
+     */
+    public double gValueDataSetsComparison(final long[] observed1, final long[] observed2)
+            throws DimensionMismatchException, NotPositiveException, ZeroException {
+
+        // Make sure lengths are same
+        if (observed1.length < 2) {
+            throw new DimensionMismatchException(observed1.length, 2);
+        }
+        if (observed1.length != observed2.length) {
+            throw new DimensionMismatchException(observed1.length, observed2.length);
+        }
+
+        // Ensure non-negative counts
+        MathArrays.checkNonNegative(observed1);
+        MathArrays.checkNonNegative(observed2);
+
+        // Compute and compare count sums
+        long countSum1 = 0;
+        long countSum2 = 0;
+
+        // Compute and compare count sums
+        final long[] collSums = new long[observed1.length];
+        final long[][] k = new long[2][observed1.length];
+
+        for (int i = 0; i < observed1.length; i++) {
+            if (observed1[i] == 0 && observed2[i] == 0) {
+                throw new ZeroException(LocalizedFormats.OBSERVED_COUNTS_BOTTH_ZERO_FOR_ENTRY,
i);
+            } else {
+                countSum1 += observed1[i];
+                countSum2 += observed2[i];
+                collSums[i] = observed1[i] + observed2[i];
+                k[0][i] = observed1[i];
+                k[1][i] = observed2[i];
+            }
+        }
+        // Ensure neither sample is uniformly 0
+        if (countSum1 == 0 || countSum2 == 0) {
+            throw new ZeroException();
+        }
+        final long[] rowSums = {countSum1, countSum2};
+        final double sum = (double) countSum1 + (double) countSum2;
+        return 2 * sum * (entropy(rowSums) + entropy(collSums) - entropy(k));
+    }
+
+    /**
+     * Calculates the root log-likelihood ratio for 2 state Datasets. See
+     * {@link #gValueDataSetsComparison(long[], long[] )}.
+     *
+     * <p>Given two events A and B, let k11 be the number of times both events
+     * occur, k12 the incidence of B without A, k21 the count of A without B,
+     * and k22 the number of times neither A nor B occurs.  What is returned
+     * by this method is </p>
+     *
+     * <p>{@code (sgn) sqrt(gValueDataSetsComparison({k11, k12}, {k21, k22})}</p>
+     *
+     * <p>where {@code sgn} is -1 if {@code k11 / (k11 + k12) < k21 / (k21 + k22))};<br/>
+     * 1 otherwise.</p>
+     *
+     * <p>Signed root LLR has two advantages over the basic LLR: a) it is positive
+     * where k11 is bigger than expected, negative where it is lower b) if there is
+     * no difference it is asymptotically normally distributed. This allows one
+     * to talk about "number of standard deviations" which is a more common frame
+     * of reference than the chi^2 distribution.</p>
+     *
+     * @param k11 number of times the two events occurred together (AB)
+     * @param k12 number of times the second event occurred WITHOUT the
+     * first event (notA,B)
+     * @param k21 number of times the first event occurred WITHOUT the
+     * second event (A, notB)
+     * @param k22 number of times something else occurred (i.e. was neither
+     * of these events (notA, notB)
+     * @return root log-likelihood ratio
+     *
+     */
+    public double rootLogLikelihoodRatio(final long k11, long k12,
+            final long k21, final long k22) {
+        final double llr = gValueDataSetsComparison(
+                new long[]{k11, k12}, new long[]{k21, k22});
+        double sqrt = FastMath.sqrt(llr);
+        if ((double) k11 / (k11 + k12) < (double) k21 / (k21 + k22)) {
+            sqrt = -sqrt;
+        }
+        return sqrt;
+    }
+
+    /**
+     * <p>Returns the <i>observed significance level</i>, or <a href=
+     * "http://www.cas.lancs.ac.uk/glossary_v1.1/hyptest.html#pvalue">
+     * p-value</a>, associated with a G-Value (Log-Likelihood Ratio) for two
+     * sample test comparing bin frequency counts in {@code observed1} and
+     * {@code observed2}.</p>
+     *
+     * <p>The number returned is the smallest significance level at which one
+     * can reject the null hypothesis that the observed counts conform to the
+     * same distribution. </p>
+     *
+     * <p>See {@link #gTestGoodnessOfFitPValue(double[], long[])} for details
+     * on how the p-value is computed.  The degrees of of freedom used to
+     * perform the test is one less than the common length of the input observed
+     * count arrays.</p>
+     *
+     * <p><strong>Preconditions</strong>:
+     * <ul> <li>Observed counts must be non-negative. </li>
+     * <li>Observed counts for a specific bin must not both be zero. </li>
+     * <li>Observed counts for a specific sample must not all be 0. </li>
+     * <li>The arrays {@code observed1} and {@ode observed2} must
+     * have the same length and their common length must be at least 2. </li>
+     * </ul><p>
+     * <p> If any of the preconditions are not met, a
+     * {@code MathIllegalArgumentException} is thrown.</p>
+     *
+     * @param observed1 array of observed frequency counts of the first data set
+     * @param observed2 array of observed frequency counts of the second data
+     * set
+     * @return p-value
+     * @throws DimensionMismatchException the the length of the arrays does not
+     * match or their common length is less than 2
+     * @throws NotPositiveException if any of the entries in {@code observed1} or
+     * {@code observed2} are negative
+     * @throws ZeroException if either all counts of {@code observed1} or
+     * {@code observed2} are zero, or if the count at some index is
+     * zero for both arrays
+     * @throws MaxCountExceededException if an error occurs computing the
+     * p-value.
+     */
+    public double gTestDataSetsComparisonPValue(final long[] observed1,
+            final long[] observed2)
+            throws DimensionMismatchException, NotPositiveException, ZeroException,
+            MaxCountExceededException {
+        final ChiSquaredDistribution distribution = new ChiSquaredDistribution(
+                (double) observed1.length - 1);
+        return 1 - distribution.cumulativeProbability(
+                gValueDataSetsComparison(observed1, observed2));
+    }
+
+    /**
+     * <p>Performs a G-Test (Log-Likelihood Ratio Test) comparing two binned
+     * data sets. The test evaluates the null hypothesis that the two lists
+     * of observed counts conform to the same frequency distribution, with
+     * significance level {@code alpha}. Returns true iff the null
+     * hypothesis can be rejected  with 100 * (1 - alpha) percent confidence.
+     * </p>
+     * <p>See {@link #gValueDataSetsComparison(long[], long[])} for details
+     * on the formula used to compute the G (LLR) statistic used in the test and
+     * {@link #gTestGoodnessOfFitPValue(double[], long[])} for information on how
+     * the observed significance level is computed. The degrees of of freedom used
+     * to perform the test is one less than the common length of the input observed
+     * count arrays. </p>
+     *
+     * <strong>Preconditions</strong>: <ul>
+     * <li>Observed counts must be non-negative. </li>
+     * <li>Observed counts for a specific bin must not both be zero. </li>
+     * <li>Observed counts for a specific sample must not all be 0. </li>
+     * <li>The arrays {@code observed1} and {@code observed2} must
+     * have the same length and their common length must be at least 2. </li>
+     * <li>{@code 0 < alpha < 0.5} </li></ul></p>
+     *
+     * <p>If any of the preconditions are not met, a
+     * {@code MathIllegalArgumentException} is thrown.</p>
+     *
+     * @param observed1 array of observed frequency counts of the first data set
+     * @param observed2 array of observed frequency counts of the second data
+     * set
+     * @param alpha significance level of the test
+     * @return true iff null hypothesis can be rejected with confidence 1 -
+     * alpha
+     * @throws DimensionMismatchException the the length of the arrays does not
+     * match
+     * @throws NotPositiveException if any of the entries in {@code observed1} or
+     * {@code observed2} are negative
+     * @throws ZeroException if either all counts of {@code observed1} or
+     * {@code observed2} are zero, or if the count at some index is
+     * zero for both arrays
+     * @throws OutOfRangeException if {@code alpha} is not in the range
+     * (0, 0.5]
+     * @throws MaxCountExceededException if an error occurs performing the test
+     */
+    public boolean gTestDataSetsComparison(
+            final long[] observed1,
+            final long[] observed2,
+            final double alpha)
+            throws DimensionMismatchException, NotPositiveException,
+            ZeroException, OutOfRangeException, MaxCountExceededException {
+
+        if (alpha <= 0 || alpha > 0.5) {
+            throw new OutOfRangeException(
+                    LocalizedFormats.OUT_OF_BOUND_SIGNIFICANCE_LEVEL, alpha, 0, 0.5);
+        }
+        return gTestDataSetsComparisonPValue(observed1, observed2) < alpha;
+    }
+}

Added: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/stat/inference/GTestTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/stat/inference/GTestTest.java?rev=1405620&view=auto
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math3/stat/inference/GTestTest.java
(added)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math3/stat/inference/GTestTest.java
Sun Nov  4 19:32:35 2012
@@ -0,0 +1,290 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.math3.stat.inference;
+
+import org.apache.commons.math3.exception.DimensionMismatchException;
+import org.apache.commons.math3.exception.NotPositiveException;
+import org.apache.commons.math3.exception.NotStrictlyPositiveException;
+import org.apache.commons.math3.exception.OutOfRangeException;
+import org.apache.commons.math3.exception.ZeroException;
+import org.junit.Assert;
+import org.junit.Test;
+
+/**
+ * Test cases for the GTest class.
+ *
+ * Data for the tests are from p64-69 in: McDonald, J.H. 2009. Handbook of
+ * Biological Statistics (2nd ed.). Sparky House Publishing, Baltimore,
+ * Maryland.
+ *
+ */
+public class GTestTest {
+
+    protected GTest testStatistic = new GTest();
+
+    @Test
+    public void testGTestGoodnesOfFit1() throws Exception {
+        final double[] exp = new double[]{
+            3d, 1d
+        };
+
+        final long[] obs = new long[]{
+            423, 133
+        };
+
+        Assert.assertEquals("G test statistic",
+                0.348721, testStatistic.gValueGoodnessOfFit(exp, obs), 1E-6);
+        final double p_gtgf = testStatistic.gTestGoodnessOfFitPValue(exp, obs);
+        Assert.assertEquals("g-Test p-value", 0.55483, p_gtgf, 1E-5);
+
+        Assert.assertFalse(testStatistic.gTestGoodnessOfFit(exp, obs, 0.05));
+    }
+
+    @Test
+    public void testGTestGoodnesOfFit2() throws Exception {
+        final double[] exp = new double[]{
+            0.54d, 0.40d, 0.05d, 0.01d
+        };
+
+        final long[] obs = new long[]{
+            70, 79, 3, 4
+        };
+        Assert.assertEquals("G test statistic",
+                13.144799, testStatistic.gValueGoodnessOfFit(exp, obs), 1E-6);
+        final double p_gtgf = testStatistic.gTestGoodnessOfFitPValue(exp, obs);
+        Assert.assertEquals("g-Test p-value", 0.004333, p_gtgf, 1E-5);
+
+        Assert.assertTrue(testStatistic.gTestGoodnessOfFit(exp, obs, 0.05));
+    }
+
+    @Test
+    public void testGTestGoodnesOfFit3() throws Exception {
+        final double[] exp = new double[]{
+            0.167d, 0.483d, 0.350d
+        };
+
+        final long[] obs = new long[]{
+            14, 21, 25
+        };
+
+        Assert.assertEquals("G test statistic",
+                4.5554, testStatistic.gValueGoodnessOfFit(exp, obs), 1E-4);
+        // Intrinisic (Hardy-Weinberg proportions) P-Value should be 0.033
+        final double p_gtgf = testStatistic.gTestGoodnessOfFitIntrinsicPValue(exp, obs);
+        Assert.assertEquals("g-Test p-value", 0.0328, p_gtgf, 1E-4);
+
+        Assert.assertFalse(testStatistic.gTestGoodnessOfFit(exp, obs, 0.05));
+    }
+
+    @Test
+    public void testGTestIndependance1() throws Exception {
+        final long[] obs1 = new long[]{
+            268, 199, 42
+        };
+
+        final long[] obs2 = new long[]{
+            807, 759, 184
+        };
+
+        final double g = testStatistic.gValueDataSetsComparison(obs1, obs2);
+
+        Assert.assertEquals("G test statistic",
+                7.3008170, g, 1E-6);
+        final double p_gti = testStatistic.gTestDataSetsComparisonPValue(obs1, obs2);
+
+        Assert.assertEquals("g-Test p-value", 0.0259805, p_gti, 1E-6);
+        Assert.assertTrue(testStatistic.gTestDataSetsComparison(obs1, obs2, 0.05));
+    }
+
+    @Test
+    public void testGTestIndependance2() throws Exception {
+        final long[] obs1 = new long[]{
+            127, 99, 264
+        };
+
+        final long[] obs2 = new long[]{
+            116, 67, 161
+        };
+
+        final double g = testStatistic.gValueDataSetsComparison(obs1, obs2);
+
+        Assert.assertEquals("G test statistic",
+                6.227288, g, 1E-6);
+        final double p_gti = testStatistic.gTestDataSetsComparisonPValue(obs1, obs2);
+
+        Assert.assertEquals("g-Test p-value", 0.04443, p_gti, 1E-5);
+        Assert.assertTrue(testStatistic.gTestDataSetsComparison(obs1, obs2, 0.05));
+    }
+
+    @Test
+    public void testGTestIndependance3() throws Exception {
+        final long[] obs1 = new long[]{
+            190, 149
+        };
+
+        final long[] obs2 = new long[]{
+            42, 49
+        };
+
+        final double g = testStatistic.gValueDataSetsComparison(obs1, obs2);
+        Assert.assertEquals("G test statistic",
+                2.8187, g, 1E-4);
+        final double p_gti = testStatistic.gTestDataSetsComparisonPValue(obs1, obs2);
+        Assert.assertEquals("g-Test p-value", 0.09317325, p_gti, 1E-6);
+
+        Assert.assertFalse(testStatistic.gTestDataSetsComparison(obs1, obs2, 0.05));
+    }
+
+    @Test
+    public void testGTestSetsComparisonBadCounts() {
+        long[] observed1 = {10, -1, 12, 10, 15};
+        long[] observed2 = {15, 10, 10, 15, 5};
+        try {
+            testStatistic.gTestDataSetsComparisonPValue(
+                    observed1, observed2);
+            Assert.fail("Expecting NotPositiveException - negative count");
+        } catch (NotPositiveException ex) {
+            // expected
+        }
+        long[] observed3 = {10, 0, 12, 10, 15};
+        long[] observed4 = {15, 0, 10, 15, 5};
+        try {
+            testStatistic.gTestDataSetsComparisonPValue(
+                    observed3, observed4);
+            Assert.fail("Expecting ZeroException - double 0's");
+        } catch (ZeroException ex) {
+            // expected
+        }
+        long[] observed5 = {10, 10, 12, 10, 15};
+        long[] observed6 = {0, 0, 0, 0, 0};
+        try {
+            testStatistic.gTestDataSetsComparisonPValue(
+                    observed5, observed6);
+            Assert.fail("Expecting ZeroException - vanishing counts");
+        } catch (ZeroException ex) {
+            // expected
+        }
+    }
+    
+    @Test
+    public void testUnmatchedArrays() {
+        final long[] observed = { 0, 1, 2, 3 };
+        final double[] expected = { 1, 1, 2 };
+        final long[] observed2 = {3, 4};
+        try {
+            testStatistic.gTestGoodnessOfFitPValue(expected, observed);
+            Assert.fail("arrays have different lengths, DimensionMismatchException expected");
+        } catch (DimensionMismatchException ex) {
+            // expected
+        }
+        try {
+            testStatistic.gTestDataSetsComparisonPValue(observed, observed2);
+            Assert.fail("arrays have different lengths, DimensionMismatchException expected");
+        } catch (DimensionMismatchException ex) {
+            // expected
+        }
+    }
+    
+    @Test
+    public void testNegativeObservedCounts() {
+        final long[] observed = { 0, 1, 2, -3 };
+        final double[] expected = { 1, 1, 2, 3};
+        final long[] observed2 = {3, 4, 5, 0};
+        try {
+            testStatistic.gTestGoodnessOfFitPValue(expected, observed);
+            Assert.fail("negative observed count, NotPositiveException expected");
+        } catch (NotPositiveException ex) {
+            // expected
+        }
+        try {
+            testStatistic.gTestDataSetsComparisonPValue(observed, observed2);
+            Assert.fail("negative observed count, NotPositiveException expected");
+        } catch (NotPositiveException ex) {
+            // expected
+        } 
+    }
+    
+    @Test
+    public void testZeroExpectedCounts() {
+        final long[] observed = { 0, 1, 2, -3 };
+        final double[] expected = { 1, 0, 2, 3};
+        try {
+            testStatistic.gTestGoodnessOfFitPValue(expected, observed);
+            Assert.fail("zero expected count, NotStrictlyPositiveException expected");
+        } catch (NotStrictlyPositiveException ex) {
+            // expected
+        }
+    }
+    
+    @Test
+    public void testBadAlpha() {
+        final long[] observed = { 0, 1, 2, 3 };
+        final double[] expected = { 1, 2, 2, 3};
+        final long[] observed2 = { 0, 2, 2, 3 };
+        try {
+            testStatistic.gTestGoodnessOfFit(expected, observed, 0.8);
+            Assert.fail("zero expected count, NotStrictlyPositiveException expected");
+        } catch (OutOfRangeException ex) {
+            // expected
+        }
+        try {
+            testStatistic.gTestDataSetsComparison(observed, observed2, -0.5);
+            Assert.fail("zero expected count, NotStrictlyPositiveException expected");
+        } catch (OutOfRangeException ex) {
+            // expected
+        }  
+    }
+    
+    @Test
+    public void testScaling() {
+      final long[] observed = {9, 11, 10, 8, 12};
+      final double[] expected1 = {10, 10, 10, 10, 10};
+      final double[] expected2 = {1000, 1000, 1000, 1000, 1000};
+      final double[] expected3 = {1, 1, 1, 1, 1};
+      final double tol = 1E-15;
+      Assert.assertEquals(
+              testStatistic.gTestGoodnessOfFitPValue(expected1, observed),
+              testStatistic.gTestGoodnessOfFitPValue(expected2, observed),
+              tol);
+      Assert.assertEquals(
+              testStatistic.gTestGoodnessOfFitPValue(expected1, observed),
+              testStatistic.gTestGoodnessOfFitPValue(expected3, observed),
+              tol);
+    }
+
+    @Test
+    public void testRootLogLikelihood() {
+        // positive where k11 is bigger than expected.
+        Assert.assertTrue(testStatistic.rootLogLikelihoodRatio(904, 21060, 1144, 283012)
> 0.0);
+
+        // negative because k11 is lower than expected
+        Assert.assertTrue(testStatistic.rootLogLikelihoodRatio(36, 21928, 60280, 623876)
< 0.0);
+
+        Assert.assertEquals(Math.sqrt(2.772589), testStatistic.rootLogLikelihoodRatio(1,
0, 0, 1), 0.000001);
+        Assert.assertEquals(-Math.sqrt(2.772589), testStatistic.rootLogLikelihoodRatio(0,
1, 1, 0), 0.000001);
+        Assert.assertEquals(Math.sqrt(27.72589), testStatistic.rootLogLikelihoodRatio(10,
0, 0, 10), 0.00001);
+
+        Assert.assertEquals(Math.sqrt(39.33052), testStatistic.rootLogLikelihoodRatio(5,
1995, 0, 100000), 0.00001);
+        Assert.assertEquals(-Math.sqrt(39.33052), testStatistic.rootLogLikelihoodRatio(0,
100000, 5, 1995), 0.00001);
+
+        Assert.assertEquals(Math.sqrt(4730.737), testStatistic.rootLogLikelihoodRatio(1000,
1995, 1000, 100000), 0.001);
+        Assert.assertEquals(-Math.sqrt(4730.737), testStatistic.rootLogLikelihoodRatio(1000,
100000, 1000, 1995), 0.001);
+
+        Assert.assertEquals(Math.sqrt(5734.343), testStatistic.rootLogLikelihoodRatio(1000,
1000, 1000, 100000), 0.001);
+        Assert.assertEquals(Math.sqrt(5714.932), testStatistic.rootLogLikelihoodRatio(1000,
1000, 1000, 99000), 0.001);
+    }
+}



Mime
View raw message