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 ≥ 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 ≥ 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 ≥ 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 * ∑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 }