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