View Javadoc
1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  package org.apache.commons.math3.stat.inference;
18  
19  import org.apache.commons.math3.distribution.ChiSquaredDistribution;
20  import org.apache.commons.math3.exception.DimensionMismatchException;
21  import org.apache.commons.math3.exception.MaxCountExceededException;
22  import org.apache.commons.math3.exception.NotPositiveException;
23  import org.apache.commons.math3.exception.NotStrictlyPositiveException;
24  import org.apache.commons.math3.exception.OutOfRangeException;
25  import org.apache.commons.math3.exception.ZeroException;
26  import org.apache.commons.math3.exception.util.LocalizedFormats;
27  import org.apache.commons.math3.util.FastMath;
28  import org.apache.commons.math3.util.MathArrays;
29  
30  /**
31   * Implements <a href="http://en.wikipedia.org/wiki/G-test">G Test</a>
32   * statistics.
33   *
34   * <p>This is known in statistical genetics as the McDonald-Kreitman test.
35   * The implementation handles both known and unknown distributions.</p>
36   *
37   * <p>Two samples tests can be used when the distribution is unknown <i>a priori</i>
38   * but provided by one sample, or when the hypothesis under test is that the two
39   * samples come from the same underlying distribution.</p>
40   *
41   * @since 3.1
42   */
43  public class GTest {
44  
45      /**
46       * Computes the <a href="http://en.wikipedia.org/wiki/G-test">G statistic
47       * for Goodness of Fit</a> comparing {@code observed} and {@code expected}
48       * frequency counts.
49       *
50       * <p>This statistic can be used to perform a G test (Log-Likelihood Ratio
51       * Test) evaluating the null hypothesis that the observed counts follow the
52       * expected distribution.</p>
53       *
54       * <p><strong>Preconditions</strong>: <ul>
55       * <li>Expected counts must all be positive. </li>
56       * <li>Observed counts must all be &ge; 0. </li>
57       * <li>The observed and expected arrays must have the same length and their
58       * common length must be at least 2. </li></ul></p>
59       *
60       * <p>If any of the preconditions are not met, a
61       * {@code MathIllegalArgumentException} is thrown.</p>
62       *
63       * <p><strong>Note:</strong>This implementation rescales the
64       * {@code expected} array if necessary to ensure that the sum of the
65       * expected and observed counts are equal.</p>
66       *
67       * @param observed array of observed frequency counts
68       * @param expected array of expected frequency counts
69       * @return G-Test statistic
70       * @throws NotPositiveException if {@code observed} has negative entries
71       * @throws NotStrictlyPositiveException if {@code expected} has entries that
72       * are not strictly positive
73       * @throws DimensionMismatchException if the array lengths do not match or
74       * are less than 2.
75       */
76      public double g(final double[] expected, final long[] observed)
77              throws NotPositiveException, NotStrictlyPositiveException,
78              DimensionMismatchException {
79  
80          if (expected.length < 2) {
81              throw new DimensionMismatchException(expected.length, 2);
82          }
83          if (expected.length != observed.length) {
84              throw new DimensionMismatchException(expected.length, observed.length);
85          }
86          MathArrays.checkPositive(expected);
87          MathArrays.checkNonNegative(observed);
88  
89          double sumExpected = 0d;
90          double sumObserved = 0d;
91          for (int i = 0; i < observed.length; i++) {
92              sumExpected += expected[i];
93              sumObserved += observed[i];
94          }
95          double ratio = 1d;
96          boolean rescale = false;
97          if (FastMath.abs(sumExpected - sumObserved) > 10E-6) {
98              ratio = sumObserved / sumExpected;
99              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</a> 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></p>
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(null, 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(null, 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)} < 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></p>
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></p>
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 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><p>
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 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(null, (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></p>
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 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 }