    001    /*
    002     * Licensed to the Apache Software Foundation (ASF) under one or more
    003     * contributor license agreements.  See the NOTICE file distributed with
    004     * this work for additional information regarding copyright ownership.
    005     * The ASF licenses this file to You under the Apache License, Version 2.0
    006     * (the "License"); you may not use this file except in compliance with
    007     * the License.  You may obtain a copy of the License at
    008     *
    009     *      http://www.apache.org/licenses/LICENSE-2.0
    010     *
    011     * Unless required by applicable law or agreed to in writing, software
    012     * distributed under the License is distributed on an "AS IS" BASIS,
    013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
    014     * See the License for the specific language governing permissions and
    015     * limitations under the License.
    016     */
    017    package org.apache.commons.math3.stat.inference;
    018
    019    import org.apache.commons.math3.distribution.ChiSquaredDistribution;
    020    import org.apache.commons.math3.exception.DimensionMismatchException;
    021    import org.apache.commons.math3.exception.MaxCountExceededException;
    022    import org.apache.commons.math3.exception.NotPositiveException;
    023    import org.apache.commons.math3.exception.NotStrictlyPositiveException;
    024    import org.apache.commons.math3.exception.OutOfRangeException;
    025    import org.apache.commons.math3.exception.ZeroException;
    026    import org.apache.commons.math3.exception.util.LocalizedFormats;
    027    import org.apache.commons.math3.util.FastMath;
    028    import org.apache.commons.math3.util.MathArrays;
    029
    030    /**
    031     * Implements <a href="http://en.wikipedia.org/wiki/G-test">G Test</a>
    032     * statistics.
    033     *
    034     * <p>This is known in statistical genetics as the McDonald-Kreitman test.
    035     * The implementation handles both known and unknown distributions.</p>
    036     *
    037     * <p>Two samples tests can be used when the distribution is unknown <i>a priori</i>
    038     * but provided by one sample, or when the hypothesis under test is that the two
    039     * samples come from the same underlying distribution.</p>
    040     *
    041     * @version $Id: GTest.java 1416643 2012-12-03 19:37:14Z tn $
    042     * @since 3.1
    043     */
    044    public class GTest {
    045
    046        /**
    047         * Computes the <a href="http://en.wikipedia.org/wiki/G-test">G statistic
    048         * for Goodness of Fit</a> comparing {@code observed} and {@code expected}
    049         * frequency counts.
    050         *
    051         * <p>This statistic can be used to perform a G test (Log-Likelihood Ratio
    052         * Test) evaluating the null hypothesis that the observed counts follow the
    053         * expected distribution.</p>
    054         *
    055         * <p><strong>Preconditions</strong>: <ul>
    056         * <li>Expected counts must all be positive. </li>
    057         * <li>Observed counts must all be &ge; 0. </li>
    058         * <li>The observed and expected arrays must have the same length and their
    059         * common length must be at least 2. </li></ul></p>
    060         *
    061         * <p>If any of the preconditions are not met, a
    062         * {@code MathIllegalArgumentException} is thrown.</p>
    063         *
    064         * <p><strong>Note:</strong>This implementation rescales the
    065         * {@code expected} array if necessary to ensure that the sum of the
    066         * expected and observed counts are equal.</p>
    067         *
    068         * @param observed array of observed frequency counts
    069         * @param expected array of expected frequency counts
    070         * @return G-Test statistic
    071         * @throws NotPositiveException if {@code observed} has negative entries
    072         * @throws NotStrictlyPositiveException if {@code expected} has entries that
    073         * are not strictly positive
    074         * @throws DimensionMismatchException if the array lengths do not match or
    075         * are less than 2.
    076         */
    077        public double g(final double[] expected, final long[] observed)
    078                throws NotPositiveException, NotStrictlyPositiveException,
    079                DimensionMismatchException {
    080
    081            if (expected.length < 2) {
    082                throw new DimensionMismatchException(expected.length, 2);
    083            }
    084            if (expected.length != observed.length) {
    085                throw new DimensionMismatchException(expected.length, observed.length);
    086            }
    087            MathArrays.checkPositive(expected);
    088            MathArrays.checkNonNegative(observed);
    089
    090            double sumExpected = 0d;
    091            double sumObserved = 0d;
    092            for (int i = 0; i < observed.length; i++) {
    093                sumExpected += expected[i];
    094                sumObserved += observed[i];
    095            }
    096            double ratio = 1d;
    097            boolean rescale = false;
    098            if (Math.abs(sumExpected - sumObserved) > 10E-6) {
    099                ratio = sumObserved / sumExpected;
    100                rescale = true;
    101            }
    102            double sum = 0d;
    103            for (int i = 0; i < observed.length; i++) {
    104                final double dev = rescale ?
    105                        FastMath.log((double) observed[i] / (ratio * expected[i])) :
    106                            FastMath.log((double) observed[i] / expected[i]);
    107                sum += ((double) observed[i]) * dev;
    108            }
    109            return 2d * sum;
    110        }
    111
    112        /**
    113         * Returns the <i>observed significance level</i>, or <a href=
    114         * "http://www.cas.lancs.ac.uk/glossary_v1.1/hyptest.html#pvalue"> p-value</a>,
    115         * associated with a G-Test for goodness of fit</a> comparing the
    116         * {@code observed} frequency counts to those in the {@code expected} array.
    117         *
    118         * <p>The number returned is the smallest significance level at which one
    119         * can reject the null hypothesis that the observed counts conform to the
    120         * frequency distribution described by the expected counts.</p>
    121         *
    122         * <p>The probability returned is the tail probability beyond
    123         * {@link #g(double[], long[]) g(expected, observed)}
    124         * in the ChiSquare distribution with degrees of freedom one less than the
    125         * common length of {@code expected} and {@code observed}.</p>
    126         *
    127         * <p> <strong>Preconditions</strong>: <ul>
    128         * <li>Expected counts must all be positive. </li>
    129         * <li>Observed counts must all be &ge; 0. </li>
    130         * <li>The observed and expected arrays must have the
    131         * same length and their common length must be at least 2.</li>
    132         * </ul></p>
    133         *
    134         * <p>If any of the preconditions are not met, a
    135         * {@code MathIllegalArgumentException} is thrown.</p>
    136         *
    137         * <p><strong>Note:</strong>This implementation rescales the
    138         * {@code expected} array if necessary to ensure that the sum of the
    139         *  expected and observed counts are equal.</p>
    140         *
    141         * @param observed array of observed frequency counts
    142         * @param expected array of expected frequency counts
    143         * @return p-value
    144         * @throws NotPositiveException if {@code observed} has negative entries
    145         * @throws NotStrictlyPositiveException if {@code expected} has entries that
    146         * are not strictly positive
    147         * @throws DimensionMismatchException if the array lengths do not match or
    148         * are less than 2.
    149         * @throws MaxCountExceededException if an error occurs computing the
    150         * p-value.
    151         */
    152        public double gTest(final double[] expected, final long[] observed)
    153                throws NotPositiveException, NotStrictlyPositiveException,
    154                DimensionMismatchException, MaxCountExceededException {
    155
    156            final ChiSquaredDistribution distribution =
    157                    new ChiSquaredDistribution(expected.length - 1.0);
    158            return 1.0 - distribution.cumulativeProbability(
    159                    g(expected, observed));
    160        }
    161
    162        /**
    163         * Returns the intrinsic (Hardy-Weinberg proportions) p-Value, as described
    164         * in p64-69 of McDonald, J.H. 2009. Handbook of Biological Statistics
    165         * (2nd ed.). Sparky House Publishing, Baltimore, Maryland.
    166         *
    167         * <p> The probability returned is the tail probability beyond
    168         * {@link #g(double[], long[]) g(expected, observed)}
    169         * in the ChiSquare distribution with degrees of freedom two less than the
    170         * common length of {@code expected} and {@code observed}.</p>
    171         *
    172         * @param observed array of observed frequency counts
    173         * @param expected array of expected frequency counts
    174         * @return p-value
    175         * @throws NotPositiveException if {@code observed} has negative entries
    176         * @throws NotStrictlyPositiveException {@code expected} has entries that are
    177         * not strictly positive
    178         * @throws DimensionMismatchException if the array lengths do not match or
    179         * are less than 2.
    180         * @throws MaxCountExceededException if an error occurs computing the
    181         * p-value.
    182         */
    183        public double gTestIntrinsic(final double[] expected, final long[] observed)
    184                throws NotPositiveException, NotStrictlyPositiveException,
    185                DimensionMismatchException, MaxCountExceededException {
    186
    187            final ChiSquaredDistribution distribution =
    188                    new ChiSquaredDistribution(expected.length - 2.0);
    189            return 1.0 - distribution.cumulativeProbability(
    190                    g(expected, observed));
    191        }
    192
    193        /**
    194         * Performs a G-Test (Log-Likelihood Ratio Test) for goodness of fit
    195         * evaluating the null hypothesis that the observed counts conform to the
    196         * frequency distribution described by the expected counts, with
    197         * significance level {@code alpha}. Returns true iff the null
    198         * hypothesis can be rejected with {@code 100 * (1 - alpha)} percent confidence.
    199         *
    200         * <p><strong>Example:</strong><br> To test the hypothesis that
    201         * {@code observed} follows {@code expected} at the 99% level,
    202         * use </p><p>
    203         * {@code gTest(expected, observed, 0.01)}</p>
    204         *
    205         * <p>Returns true iff {@link #gTest(double[], long[])
    206         *  gTestGoodnessOfFitPValue(expected, observed)} < alpha</p>
    207         *
    208         * <p><strong>Preconditions</strong>: <ul>
    209         * <li>Expected counts must all be positive. </li>
    210         * <li>Observed counts must all be &ge; 0. </li>
    211         * <li>The observed and expected arrays must have the same length and their
    212         * common length must be at least 2.
    213         * <li> {@code 0 < alpha < 0.5} </li></ul></p>
    214         *
    215         * <p>If any of the preconditions are not met, a
    216         * {@code MathIllegalArgumentException} is thrown.</p>
    217         *
    218         * <p><strong>Note:</strong>This implementation rescales the
    219         * {@code expected} array if necessary to ensure that the sum of the
    220         * expected and observed counts are equal.</p>
    221         *
    222         * @param observed array of observed frequency counts
    223         * @param expected array of expected frequency counts
    224         * @param alpha significance level of the test
    225         * @return true iff null hypothesis can be rejected with confidence 1 -
    226         * alpha
    227         * @throws NotPositiveException if {@code observed} has negative entries
    228         * @throws NotStrictlyPositiveException if {@code expected} has entries that
    229         * are not strictly positive
    230         * @throws DimensionMismatchException if the array lengths do not match or
    231         * are less than 2.
    232         * @throws MaxCountExceededException if an error occurs computing the
    233         * p-value.
    234         * @throws OutOfRangeException if alpha is not strictly greater than zero
    235         * and less than or equal to 0.5
    236         */
    237        public boolean gTest(final double[] expected, final long[] observed,
    238                final double alpha)
    239                throws NotPositiveException, NotStrictlyPositiveException,
    240                DimensionMismatchException, OutOfRangeException, MaxCountExceededException {
    241
    242            if ((alpha <= 0) || (alpha > 0.5)) {
    243                throw new OutOfRangeException(LocalizedFormats.OUT_OF_BOUND_SIGNIFICANCE_LEVEL,
    244                        alpha, 0, 0.5);
    245            }
    246            return gTest(expected, observed) < alpha;
    247        }
    248
    249        /**
    250         * Calculates the <a href=
    251         * "http://en.wikipedia.org/wiki/Entropy_%28information_theory%29">Shannon
    252         * entropy</a> for 2 Dimensional Matrix.  The value returned is the entropy
    253         * of the vector formed by concatenating the rows (or columns) of {@code k}
    254         * to form a vector. See {@link #entropy(long[])}.
    255         *
    256         * @param k 2 Dimensional Matrix of long values (for ex. the counts of a
    257         * trials)
    258         * @return Shannon Entropy of the given Matrix
    259         *
    260         */
    261        private double entropy(final long[][] k) {
    262            double h = 0d;
    263            double sum_k = 0d;
    264            for (int i = 0; i < k.length; i++) {
    265                for (int j = 0; j < k[i].length; j++) {
    266                    sum_k += (double) k[i][j];
    267                }
    268            }
    269            for (int i = 0; i < k.length; i++) {
    270                for (int j = 0; j < k[i].length; j++) {
    271                    if (k[i][j] != 0) {
    272                        final double p_ij = (double) k[i][j] / sum_k;
    273                        h += p_ij * Math.log(p_ij);
    274                    }
    275                }
    276            }
    277            return -h;
    278        }
    279
    280        /**
    281         * Calculates the <a href="http://en.wikipedia.org/wiki/Entropy_%28information_theory%29">
    282         * Shannon entropy</a> for a vector.  The values of {@code k} are taken to be
    283         * incidence counts of the values of a random variable. What is returned is <br/>
    284         * &sum;p<sub>i</sub>log(p<sub>i</sub><br/>
    285         * where p<sub>i</sub> = k[i] / (sum of elements in k)
    286         *
    287         * @param k Vector (for ex. Row Sums of a trials)
    288         * @return Shannon Entropy of the given Vector
    289         *
    290         */
    291        private double entropy(final long[] k) {
    292            double h = 0d;
    293            double sum_k = 0d;
    294            for (int i = 0; i < k.length; i++) {
    295                sum_k += (double) k[i];
    296            }
    297            for (int i = 0; i < k.length; i++) {
    298                if (k[i] != 0) {
    299                    final double p_i = (double) k[i] / sum_k;
    300                    h += p_i * Math.log(p_i);
    301                }
    302            }
    303            return -h;
    304        }
    305
    306        /**
    307         * <p>Computes a G (Log-Likelihood Ratio) two sample test statistic for
    308         * independence comparing frequency counts in
    309         * {@code observed1} and {@code observed2}. The sums of frequency
    310         * counts in the two samples are not required to be the same. The formula
    311         * used to compute the test statistic is </p>
    312         *
    313         * <p>{@code 2 * totalSum * [H(rowSums) + H(colSums) - H(k)]}</p>
    314         *
    315         * <p> where {@code H} is the
    316         * <a href="http://en.wikipedia.org/wiki/Entropy_%28information_theory%29">
    317         * Shannon Entropy</a> of the random variable formed by viewing the elements
    318         * of the argument array as incidence counts; <br/>
    319         * {@code k} is a matrix with rows {@code [observed1, observed2]}; <br/>
    320         * {@code rowSums, colSums} are the row/col sums of {@code k}; <br>
    321         * and {@code totalSum} is the overall sum of all entries in {@code k}.</p>
    322         *
    323         * <p>This statistic can be used to perform a G test evaluating the null
    324         * hypothesis that both observed counts are independent </p>
    325         *
    326         * <p> <strong>Preconditions</strong>: <ul>
    327         * <li>Observed counts must be non-negative. </li>
    328         * <li>Observed counts for a specific bin must not both be zero. </li>
    329         * <li>Observed counts for a specific sample must not all be  0. </li>
    330         * <li>The arrays {@code observed1} and {@code observed2} must have
    331         * the same length and their common length must be at least 2. </li></ul></p>
    332         *
    333         * <p>If any of the preconditions are not met, a
    334         * {@code MathIllegalArgumentException} is thrown.</p>
    335         *
    336         * @param observed1 array of observed frequency counts of the first data set
    337         * @param observed2 array of observed frequency counts of the second data
    338         * set
    339         * @return G-Test statistic
    340         * @throws DimensionMismatchException the the lengths of the arrays do not
    341         * match or their common length is less than 2
    342         * @throws NotPositiveException if any entry in {@code observed1} or
    343         * {@code observed2} is negative
    344         * @throws ZeroException if either all counts of
    345         * {@code observed1} or {@code observed2} are zero, or if the count
    346         * at the same index is zero for both arrays.
    347         */
    348        public double gDataSetsComparison(final long[] observed1, final long[] observed2)
    349                throws DimensionMismatchException, NotPositiveException, ZeroException {
    350
    351            // Make sure lengths are same
    352            if (observed1.length < 2) {
    353                throw new DimensionMismatchException(observed1.length, 2);
    354            }
    355            if (observed1.length != observed2.length) {
    356                throw new DimensionMismatchException(observed1.length, observed2.length);
    357            }
    358
    359            // Ensure non-negative counts
    360            MathArrays.checkNonNegative(observed1);
    361            MathArrays.checkNonNegative(observed2);
    362
    363            // Compute and compare count sums
    364            long countSum1 = 0;
    365            long countSum2 = 0;
    366
    367            // Compute and compare count sums
    368            final long[] collSums = new long[observed1.length];
    369            final long[][] k = new long[2][observed1.length];
    370
    371            for (int i = 0; i < observed1.length; i++) {
    372                if (observed1[i] == 0 && observed2[i] == 0) {
    373                    throw new ZeroException(LocalizedFormats.OBSERVED_COUNTS_BOTTH_ZERO_FOR_ENTRY, i);
    374                } else {
    375                    countSum1 += observed1[i];
    376                    countSum2 += observed2[i];
    377                    collSums[i] = observed1[i] + observed2[i];
    378                    k[0][i] = observed1[i];
    379                    k[1][i] = observed2[i];
    380                }
    381            }
    382            // Ensure neither sample is uniformly 0
    383            if (countSum1 == 0 || countSum2 == 0) {
    384                throw new ZeroException();
    385            }
    386            final long[] rowSums = {countSum1, countSum2};
    387            final double sum = (double) countSum1 + (double) countSum2;
    388            return 2 * sum * (entropy(rowSums) + entropy(collSums) - entropy(k));
    389        }
    390
    391        /**
    392         * Calculates the root log-likelihood ratio for 2 state Datasets. See
    393         * {@link #gDataSetsComparison(long[], long[] )}.
    394         *
    395         * <p>Given two events A and B, let k11 be the number of times both events
    396         * occur, k12 the incidence of B without A, k21 the count of A without B,
    397         * and k22 the number of times neither A nor B occurs.  What is returned
    398         * by this method is </p>
    399         *
    400         * <p>{@code (sgn) sqrt(gValueDataSetsComparison({k11, k12}, {k21, k22})}</p>
    401         *
    402         * <p>where {@code sgn} is -1 if {@code k11 / (k11 + k12) < k21 / (k21 + k22))};<br/>
    403         * 1 otherwise.</p>
    404         *
    405         * <p>Signed root LLR has two advantages over the basic LLR: a) it is positive
    406         * where k11 is bigger than expected, negative where it is lower b) if there is
    407         * no difference it is asymptotically normally distributed. This allows one
    408         * to talk about "number of standard deviations" which is a more common frame
    409         * of reference than the chi^2 distribution.</p>
    410         *
    411         * @param k11 number of times the two events occurred together (AB)
    412         * @param k12 number of times the second event occurred WITHOUT the
    413         * first event (notA,B)
    414         * @param k21 number of times the first event occurred WITHOUT the
    415         * second event (A, notB)
    416         * @param k22 number of times something else occurred (i.e. was neither
    417         * of these events (notA, notB)
    418         * @return root log-likelihood ratio
    419         *
    420         */
    421        public double rootLogLikelihoodRatio(final long k11, long k12,
    422                final long k21, final long k22) {
    423            final double llr = gDataSetsComparison(
    424                    new long[]{k11, k12}, new long[]{k21, k22});
    425            double sqrt = FastMath.sqrt(llr);
    426            if ((double) k11 / (k11 + k12) < (double) k21 / (k21 + k22)) {
    427                sqrt = -sqrt;
    428            }
    429            return sqrt;
    430        }
    431
    432        /**
    433         * <p>Returns the <i>observed significance level</i>, or <a href=
    434         * "http://www.cas.lancs.ac.uk/glossary_v1.1/hyptest.html#pvalue">
    435         * p-value</a>, associated with a G-Value (Log-Likelihood Ratio) for two
    436         * sample test comparing bin frequency counts in {@code observed1} and
    437         * {@code observed2}.</p>
    438         *
    439         * <p>The number returned is the smallest significance level at which one
    440         * can reject the null hypothesis that the observed counts conform to the
    441         * same distribution. </p>
    442         *
    443         * <p>See {@link #gTest(double[], long[])} for details
    444         * on how the p-value is computed.  The degrees of of freedom used to
    445         * perform the test is one less than the common length of the input observed
    446         * count arrays.</p>
    447         *
    448         * <p><strong>Preconditions</strong>:
    449         * <ul> <li>Observed counts must be non-negative. </li>
    450         * <li>Observed counts for a specific bin must not both be zero. </li>
    451         * <li>Observed counts for a specific sample must not all be 0. </li>
    452         * <li>The arrays {@code observed1} and {@code observed2} must
    453         * have the same length and their common length must be at least 2. </li>
    454         * </ul><p>
    455         * <p> If any of the preconditions are not met, a
    456         * {@code MathIllegalArgumentException} is thrown.</p>
    457         *
    458         * @param observed1 array of observed frequency counts of the first data set
    459         * @param observed2 array of observed frequency counts of the second data
    460         * set
    461         * @return p-value
    462         * @throws DimensionMismatchException the the length of the arrays does not
    463         * match or their common length is less than 2
    464         * @throws NotPositiveException if any of the entries in {@code observed1} or
    465         * {@code observed2} are negative
    466         * @throws ZeroException if either all counts of {@code observed1} or
    467         * {@code observed2} are zero, or if the count at some index is
    468         * zero for both arrays
    469         * @throws MaxCountExceededException if an error occurs computing the
    470         * p-value.
    471         */
    472        public double gTestDataSetsComparison(final long[] observed1,
    473                final long[] observed2)
    474                throws DimensionMismatchException, NotPositiveException, ZeroException,
    475                MaxCountExceededException {
    476            final ChiSquaredDistribution distribution = new ChiSquaredDistribution(
    477                    (double) observed1.length - 1);
    478            return 1 - distribution.cumulativeProbability(
    479                    gDataSetsComparison(observed1, observed2));
    480        }
    481
    482        /**
    483         * <p>Performs a G-Test (Log-Likelihood Ratio Test) comparing two binned
    484         * data sets. The test evaluates the null hypothesis that the two lists
    485         * of observed counts conform to the same frequency distribution, with
    486         * significance level {@code alpha}. Returns true iff the null
    487         * hypothesis can be rejected  with 100 * (1 - alpha) percent confidence.
    488         * </p>
    489         * <p>See {@link #gDataSetsComparison(long[], long[])} for details
    490         * on the formula used to compute the G (LLR) statistic used in the test and
    491         * {@link #gTest(double[], long[])} for information on how
    492         * the observed significance level is computed. The degrees of of freedom used
    493         * to perform the test is one less than the common length of the input observed
    494         * count arrays. </p>
    495         *
    496         * <strong>Preconditions</strong>: <ul>
    497         * <li>Observed counts must be non-negative. </li>
    498         * <li>Observed counts for a specific bin must not both be zero. </li>
    499         * <li>Observed counts for a specific sample must not all be 0. </li>
    500         * <li>The arrays {@code observed1} and {@code observed2} must
    501         * have the same length and their common length must be at least 2. </li>
    502         * <li>{@code 0 < alpha < 0.5} </li></ul></p>
    503         *
    504         * <p>If any of the preconditions are not met, a
    505         * {@code MathIllegalArgumentException} is thrown.</p>
    506         *
    507         * @param observed1 array of observed frequency counts of the first data set
    508         * @param observed2 array of observed frequency counts of the second data
    509         * set
    510         * @param alpha significance level of the test
    511         * @return true iff null hypothesis can be rejected with confidence 1 -
    512         * alpha
    513         * @throws DimensionMismatchException the the length of the arrays does not
    514         * match
    515         * @throws NotPositiveException if any of the entries in {@code observed1} or
    516         * {@code observed2} are negative
    517         * @throws ZeroException if either all counts of {@code observed1} or
    518         * {@code observed2} are zero, or if the count at some index is
    519         * zero for both arrays
    520         * @throws OutOfRangeException if {@code alpha} is not in the range
    521         * (0, 0.5]
    522         * @throws MaxCountExceededException if an error occurs performing the test
    523         */
    524        public boolean gTestDataSetsComparison(
    525                final long[] observed1,
    526                final long[] observed2,
    527                final double alpha)
    528                throws DimensionMismatchException, NotPositiveException,
    529                ZeroException, OutOfRangeException, MaxCountExceededException {
    530
    531            if (alpha <= 0 || alpha > 0.5) {
    532                throw new OutOfRangeException(
    533                        LocalizedFormats.OUT_OF_BOUND_SIGNIFICANCE_LEVEL, alpha, 0, 0.5);
    534            }
    535            return gTestDataSetsComparison(observed1, observed2) < alpha;
    536        }
    537    }