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.math3.stat.inference; 018 019import org.apache.commons.math3.distribution.ChiSquaredDistribution; 020import org.apache.commons.math3.exception.DimensionMismatchException; 021import org.apache.commons.math3.exception.MaxCountExceededException; 022import org.apache.commons.math3.exception.NotPositiveException; 023import org.apache.commons.math3.exception.NotStrictlyPositiveException; 024import org.apache.commons.math3.exception.OutOfRangeException; 025import org.apache.commons.math3.exception.ZeroException; 026import org.apache.commons.math3.exception.util.LocalizedFormats; 027import org.apache.commons.math3.util.FastMath; 028import 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 * @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 ≥ 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></p> 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</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 ≥ 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 ≥ 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 * ∑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}