ContingencyTables.java
Upload User: rhdiban
Upload Date: 2013-08-09
Package Size: 15085k
Code Size: 18k
Category:

Windows Develop

Development Platform:

Java

  1. /*
  2.  *    This program is free software; you can redistribute it and/or modify
  3.  *    it under the terms of the GNU General Public License as published by
  4.  *    the Free Software Foundation; either version 2 of the License, or
  5.  *    (at your option) any later version.
  6.  *
  7.  *    This program is distributed in the hope that it will be useful,
  8.  *    but WITHOUT ANY WARRANTY; without even the implied warranty of
  9.  *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  10.  *    GNU General Public License for more details.
  11.  *
  12.  *    You should have received a copy of the GNU General Public License
  13.  *    along with this program; if not, write to the Free Software
  14.  *    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  15.  */
  16. /*
  17.  *    ContingencyTables.java
  18.  *    Copyright (C) 1999 Eibe Frank
  19.  *
  20.  */
  21. package weka.core;
  22. /**
  23.  * Class implementing some statistical routines for contingency tables.
  24.  *
  25.  * @author Eibe Frank (eibe@cs.waikato.ac.nz)
  26.  * @version $Revision: 1.4 $
  27.  */
  28. public class ContingencyTables {
  29.   /** The natural logarithm of 2 */
  30.   private static double log2 = Math.log(2);
  31.   /**
  32.    * Returns chi-squared probability for a given matrix.
  33.    *
  34.    * @param matrix the contigency table
  35.    * @param yates is Yates' correction to be used?
  36.    * @return the chi-squared probability
  37.    */
  38.   public static double chiSquared(double [][] matrix, boolean yates) {
  39.     int df = (matrix.length - 1) * (matrix[0].length - 1);
  40.     return Statistics.chiSquaredProbability(chiVal(matrix, yates), df);
  41.   }
  42.   /**
  43.    * Computes chi-squared statistic for a contingency table.
  44.    *
  45.    * @param matrix the contigency table
  46.    * @param yates is Yates' correction to be used?
  47.    * @return the value of the chi-squared statistic
  48.    */
  49.   public static double chiVal(double [][] matrix, boolean useYates) {
  50.     
  51.     int df, nrows, ncols, row, col;
  52.     double[] rtotal, ctotal;
  53.     double expect = 0, chival = 0, n = 0;
  54.     boolean yates = true;
  55.     
  56.     nrows = matrix.length;
  57.     ncols = matrix[0].length;
  58.     rtotal = new double [nrows];
  59.     ctotal = new double [ncols];
  60.     for (row = 0; row < nrows; row++) {
  61.       for (col = 0; col < ncols; col++) {
  62. rtotal[row] += matrix[row][col];
  63. ctotal[col] += matrix[row][col];
  64. n += matrix[row][col];
  65.       }
  66.     }
  67.     df = (nrows - 1)*(ncols - 1);
  68.     if ((df > 1) || (!useYates)) {
  69.       yates = false;
  70.     } else if (df <= 0) {
  71.       return 0;
  72.     }
  73.     chival = 0.0;
  74.     for (row = 0; row < nrows; row++) {
  75.       if (Utils.gr(rtotal[row], 0)) {
  76. for (col = 0; col < ncols; col++) {
  77.   if (Utils.gr(ctotal[col], 0)) {
  78.     expect = (ctotal[col] * rtotal[row]) / n;
  79.     chival += chiCell (matrix[row][col], expect, yates);
  80.   }
  81. }
  82.       }
  83.     }
  84.     return chival;
  85.   }
  86.   /**
  87.    * Tests if Cochran's criterion is fullfilled for the given
  88.    * contingency table. Rows and columns with all zeros are not considered
  89.    * relevant.
  90.    *
  91.    * @param matrix the contigency table to be tested
  92.    * @return true if contingency table is ok, false if not
  93.    */
  94.   public static boolean cochransCriterion(double[][] matrix) {
  95.     double[] rtotal, ctotal;
  96.     double n = 0, expect, smallfreq = 5;
  97.     int smallcount = 0, nonZeroRows = 0, nonZeroColumns = 0, nrows, ncols, 
  98.       row, col;
  99.     nrows = matrix.length;
  100.     ncols = matrix[0].length;
  101.     rtotal = new double [nrows];
  102.     ctotal = new double [ncols];
  103.     for (row = 0; row < nrows; row++) {
  104.       for (col = 0; col < ncols; col++) {
  105. rtotal[row] += matrix[row][col];
  106. ctotal[col] += matrix[row][col];
  107. n += matrix[row][col];
  108.       }
  109.     }
  110.     for (row = 0; row < nrows; row++) {
  111.       if (Utils.gr(rtotal[row], 0)) {
  112. nonZeroRows++;
  113.       }
  114.     }
  115.     for (col = 0; col < ncols; col++) {
  116.       if (Utils.gr(ctotal[col], 0)) {
  117. nonZeroColumns++;
  118.       }
  119.     }
  120.     for (row = 0; row < nrows; row++) {
  121.       if (Utils.gr(rtotal[row], 0)) {
  122. for (col = 0; col < ncols; col++) {
  123.   if (Utils.gr(ctotal[col], 0)) {
  124.     expect = (ctotal[col] * rtotal[row]) / n;
  125.     if (Utils.sm(expect, smallfreq)) {
  126.       if (Utils.sm(expect, 1)) {
  127. return false;
  128.       } else {
  129. smallcount++;
  130. if (smallcount > (nonZeroRows * nonZeroColumns) / smallfreq) {
  131.   return false;
  132. }
  133.       }
  134.     }
  135.   }
  136. }
  137.       }
  138.     }
  139.     return true;
  140.   }
  141.   /**
  142.    * Computes Cramer's V for a contingency table.
  143.    *
  144.    * @param matrix the contingency table
  145.    * @return Cramer's V
  146.    */
  147.   public static double CramersV(double [][] matrix) {
  148.     int row, col, nrows,ncols, min;
  149.     double n = 0;
  150.     
  151.     nrows = matrix.length;
  152.     ncols = matrix[0].length;
  153.     for (row = 0; row < nrows; row++) {
  154.       for (col = 0; col < ncols; col++) {
  155. n += matrix[row][col];
  156.       }
  157.     }
  158.     min = nrows < ncols ? nrows-1 : ncols-1;
  159.     if ((min == 0) || Utils.eq(n, 0))
  160.       return 0;
  161.     return Math.sqrt(chiVal(matrix, false) / (n * (double)min)); 
  162.   } 
  163.   /**
  164.    * Computes the entropy of the given array.
  165.    *
  166.    * @param array the array
  167.    * @return the entropy
  168.    */
  169.   public static double entropy(double[] array) {
  170.     double returnValue = 0, sum = 0;
  171.     for (int i = 0; i < array.length; i++) {
  172.       returnValue -= lnFunc(array[i]);
  173.       sum += array[i];
  174.     }
  175.     if (Utils.eq(sum, 0)) {
  176.       return 0;
  177.     } else {
  178.       return (returnValue + lnFunc(sum)) / (sum * log2);
  179.     }
  180.   }
  181.   /**
  182.    * Computes conditional entropy of the rows given
  183.    * the columns.
  184.    *
  185.    * @param matrix the contingency table
  186.    * @return the conditional entropy of the rows given the columns
  187.    */
  188.   public static double entropyConditionedOnColumns(double[][] matrix) {
  189.     
  190.     double returnValue = 0, sumForColumn, total = 0;
  191.     for (int j = 0; j < matrix[0].length; j++) {
  192.       sumForColumn = 0;
  193.       for (int i = 0; i < matrix.length; i++) {
  194. returnValue = returnValue + lnFunc(matrix[i][j]);
  195. sumForColumn += matrix[i][j];
  196.       }
  197.       returnValue = returnValue - lnFunc(sumForColumn);
  198.       total += sumForColumn;
  199.     }
  200.     if (Utils.eq(total, 0)) {
  201.       return 0;
  202.     }
  203.     return -returnValue / (total * log2);
  204.   }
  205.   /**
  206.    * Computes conditional entropy of the columns given
  207.    * the rows.
  208.    *
  209.    * @param matrix the contingency table
  210.    * @return the conditional entropy of the columns given the rows
  211.    */
  212.   public static double entropyConditionedOnRows(double[][] matrix) {
  213.     
  214.     double returnValue = 0, sumForRow, total = 0;
  215.     for (int i = 0; i < matrix.length; i++) {
  216.       sumForRow = 0;
  217.       for (int j = 0; j < matrix[0].length; j++) {
  218. returnValue = returnValue + lnFunc(matrix[i][j]);
  219. sumForRow += matrix[i][j];
  220.       }
  221.       returnValue = returnValue - lnFunc(sumForRow);
  222.       total += sumForRow;
  223.     }
  224.     if (Utils.eq(total, 0)) {
  225.       return 0;
  226.     }
  227.     return -returnValue / (total * log2);
  228.   }
  229.   /**
  230.    * Computes conditional entropy of the columns given the rows
  231.    * of the test matrix with respect to the train matrix. Uses a
  232.    * Laplace prior. Does NOT normalize the entropy.
  233.    *
  234.    * @param train the train matrix 
  235.    * @param test the test matrix
  236.    * @param the number of symbols for Laplace
  237.    * @return the entropy
  238.    */
  239.   public static double entropyConditionedOnRows(double[][] train, 
  240. double[][] test,
  241. double numClasses) {
  242.     
  243.     double returnValue = 0, trainSumForRow, testSumForRow, testSum = 0;
  244.     for (int i = 0; i < test.length; i++) {
  245.       trainSumForRow = 0;
  246.       testSumForRow = 0;
  247.       for (int j = 0; j < test[0].length; j++) {
  248. returnValue -= test[i][j] * Math.log(train[i][j] + 1);
  249. trainSumForRow += train[i][j];
  250. testSumForRow += test[i][j];
  251.       }
  252.       testSum = testSumForRow;
  253.       returnValue += testSumForRow * Math.log(trainSumForRow + 
  254.      numClasses);
  255.     }
  256.     return returnValue / (testSum * log2);
  257.   }
  258.   /**
  259.    * Computes the rows' entropy for the given contingency table.
  260.    *
  261.    * @param matrix the contingency table
  262.    * @return the rows' entropy
  263.    */
  264.   public static double entropyOverRows(double[][] matrix) {
  265.     
  266.     double returnValue = 0, sumForRow, total = 0;
  267.     for (int i = 0; i < matrix.length; i++) {
  268.       sumForRow = 0;
  269.       for (int j = 0; j < matrix[0].length; j++) {
  270. sumForRow += matrix[i][j];
  271.       }
  272.       returnValue = returnValue - lnFunc(sumForRow);
  273.       total += sumForRow;
  274.     }
  275.     if (Utils.eq(total, 0)) {
  276.       return 0;
  277.     }
  278.     return (returnValue + lnFunc(total)) / (total * log2);
  279.   }
  280.   /**
  281.    * Computes the columns' entropy for the given contingency table.
  282.    *
  283.    * @param matrix the contingency table
  284.    * @return the columns' entropy
  285.    */
  286.   public static double entropyOverColumns(double[][] matrix){
  287.     
  288.     double returnValue = 0, sumForColumn, total = 0;
  289.     for (int j = 0; j < matrix[0].length; j++){
  290.       sumForColumn = 0;
  291.       for (int i = 0; i < matrix.length; i++) {
  292. sumForColumn += matrix[i][j];
  293.       }
  294.       returnValue = returnValue - lnFunc(sumForColumn);
  295.       total += sumForColumn;
  296.     }
  297.     if (Utils.eq(total, 0)) {
  298.       return 0;
  299.     }
  300.     return (returnValue + lnFunc(total)) / (total * log2);
  301.   }
  302.   /**
  303.    * Computes gain ratio for contingency table (split on rows).
  304.    * Returns Double.MAX_VALUE if the split entropy is 0.
  305.    *
  306.    * @param matrix the contingency table
  307.    * @return the gain ratio
  308.    */
  309.   public static double gainRatio(double[][] matrix){
  310.     
  311.     double preSplit = 0, postSplit = 0, splitEnt = 0,
  312.       sumForRow, sumForColumn, total = 0, infoGain;
  313.     // Compute entropy before split
  314.     for (int i = 0; i < matrix[0].length; i++) {
  315.       sumForColumn = 0;
  316.       for (int j = 0; j < matrix.length; j++) 
  317. sumForColumn += matrix[j][i];
  318.       preSplit += lnFunc(sumForColumn);
  319.       total += sumForColumn;
  320.     }
  321.     preSplit -= lnFunc(total);
  322.     // Compute entropy after split and split entropy
  323.     for (int i = 0; i < matrix.length; i++) {
  324.       sumForRow = 0;
  325.       for (int j = 0; j < matrix[0].length; j++) {
  326. postSplit += lnFunc(matrix[i][j]);
  327. sumForRow += matrix[i][j];
  328.       }
  329.       splitEnt += lnFunc(sumForRow);
  330.     }
  331.     postSplit -= splitEnt;
  332.     splitEnt -= lnFunc(total);
  333.     infoGain = preSplit - postSplit;
  334.     if (Utils.eq(splitEnt, 0))
  335.       return 0;
  336.     return infoGain / splitEnt;
  337.   }
  338.   /**
  339.    * Returns negative base 2 logarithm of multiple hypergeometric
  340.    * probability for a contingency table.
  341.    *
  342.    * @param matrix the contingency table
  343.    * @return the log of the hypergeometric probability of the contingency table 
  344.    */
  345.   public static double log2MultipleHypergeometric(double[][] matrix) {
  346.     double sum = 0, sumForRow, sumForColumn, total = 0;
  347.     for (int i = 0; i < matrix.length; i++) {
  348.       sumForRow = 0;
  349.       for (int j = 0; j < matrix[i].length; j++) {
  350. sumForRow += matrix[i][j];
  351.       }
  352.       sum += SpecialFunctions.lnFactorial(sumForRow);
  353.       total += sumForRow;
  354.     }
  355.     for (int j = 0; j < matrix[0].length; j++) {
  356.       sumForColumn = 0;
  357.       for (int i = 0; i < matrix.length; i++) {
  358. sumForColumn += matrix [i][j];
  359.       }
  360.       sum += SpecialFunctions.lnFactorial(sumForColumn);
  361.     }
  362.     for (int i = 0; i < matrix.length; i++) {
  363.       for (int j = 0; j < matrix[i].length; j++) {
  364. sum -= SpecialFunctions.lnFactorial(matrix[i][j]);
  365.       }
  366.     }
  367.     sum -= SpecialFunctions.lnFactorial(total);
  368.     return -sum / log2;
  369.   }
  370.   /**
  371.    * Reduces a matrix by deleting all zero rows and columns.
  372.    *
  373.    * @param matrix the matrix to be reduced
  374.    * @param the matrix with all zero rows and columns deleted
  375.    */
  376.   public static double[][] reduceMatrix(double[][] matrix) {
  377.     int row, col, currCol, currRow, nrows, ncols, 
  378.       nonZeroRows = 0, nonZeroColumns = 0;
  379.     double[] rtotal, ctotal;
  380.     double[][] newMatrix;
  381.     nrows = matrix.length;
  382.     ncols = matrix[0].length;
  383.     rtotal = new double [nrows];
  384.     ctotal = new double [ncols];
  385.     for (row = 0; row < nrows; row++) {
  386.       for (col = 0; col < ncols; col++) {
  387. rtotal[row] += matrix[row][col];
  388. ctotal[col] += matrix[row][col];
  389.       }
  390.     }
  391.     for (row = 0; row < nrows; row++) {
  392.       if (Utils.gr(rtotal[row],0)) {
  393. nonZeroRows++;
  394.       }
  395.     }
  396.     for (col = 0; col < ncols; col++) {
  397.       if (Utils.gr(ctotal[col],0)) {
  398. nonZeroColumns++;
  399.       }
  400.     }
  401.     newMatrix = new double[nonZeroRows][nonZeroColumns];
  402.     currRow = 0;
  403.     for (row = 0; row < nrows; row++) {
  404.       if (Utils.gr(rtotal[row],0)) {
  405. currCol = 0;
  406. for (col = 0; col < ncols; col++) {
  407.   if (Utils.gr(ctotal[col],0)) {
  408.     newMatrix[currRow][currCol] = matrix[row][col];
  409.     currCol++;
  410.   }
  411. }
  412. currRow++;
  413.       }
  414.     }
  415.     return newMatrix;
  416.   }
  417.     
  418.    /**
  419.     * Calculates the symmetrical uncertainty for base 2.
  420.     *
  421.     * @param matrix the contingency table
  422.     * @return the calculated symmetrical uncertainty
  423.     *
  424.     */
  425.   public static double symmetricalUncertainty(double matrix[][]) {
  426.     double sumForColumn, sumForRow, total = 0, columnEntropy = 0, 
  427.       rowEntropy = 0, entropyConditionedOnRows = 0, infoGain = 0;
  428.     // Compute entropy for columns
  429.     for (int i = 0; i < matrix[0].length; i++) {
  430.       sumForColumn = 0;
  431.       for (int j = 0; j < matrix.length; j++) {
  432. sumForColumn += matrix[j][i];
  433.       }
  434.       columnEntropy += lnFunc(sumForColumn);
  435.       total += sumForColumn;
  436.     }
  437.     columnEntropy -= lnFunc(total);
  438.     // Compute entropy for rows and conditional entropy
  439.     for (int i = 0; i < matrix.length; i++) {
  440.       sumForRow = 0;
  441.       for (int j = 0; j < matrix[0].length; j++) { 
  442. sumForRow += matrix[i][j];
  443. entropyConditionedOnRows += lnFunc(matrix[i][j]);
  444.       }
  445.       rowEntropy += lnFunc(sumForRow);
  446.     }
  447.     entropyConditionedOnRows -= rowEntropy;
  448.     rowEntropy -= lnFunc(total);
  449.     infoGain = columnEntropy - entropyConditionedOnRows;
  450.     if (Utils.eq(columnEntropy, 0) || Utils.eq(rowEntropy, 0))
  451.       return 0;
  452.     return 2.0 * (infoGain / (columnEntropy + rowEntropy));
  453.   }
  454.   /**
  455.    * Computes Goodman and Kruskal's tau-value for a contingency table.
  456.    *
  457.    * @param matrix the contingency table
  458.    * @param Goodman and Kruskal's tau-value
  459.    */
  460.   public static double tauVal(double[][] matrix) {
  461.     
  462.     int nrows, ncols, row, col;
  463.     double [] ctotal;
  464.     double maxcol = 0, max, maxtotal = 0, n = 0;
  465.     
  466.     nrows = matrix.length;
  467.     ncols = matrix[0].length;
  468.     ctotal = new double [ncols];
  469.     for (row = 0; row < nrows; row++) {
  470.       max = 0;
  471.       for (col = 0; col < ncols; col++) {
  472. if (Utils.gr(matrix[row][col], max)) 
  473.   max = matrix[row][col];
  474. ctotal[col] += matrix[row][col];
  475. n += matrix[row][col];
  476.       }
  477.       maxtotal += max;
  478.     }
  479.     if (Utils.eq(n, 0)) {
  480.       return 0;
  481.     }
  482.     maxcol = ctotal[Utils.maxIndex(ctotal)];
  483.     return (maxtotal - maxcol)/(n - maxcol);
  484.   }
  485.   /**
  486.    * Help method for computing entropy.
  487.    */
  488.   private static double lnFunc(double num){
  489.     
  490.     // Constant hard coded for efficiency reasons
  491.     if (num < 1e-6) {
  492.       return 0;
  493.     } else {
  494.       return num * Math.log(num);
  495.     }
  496.   }
  497.   
  498.   /**
  499.    * Computes chi-value for one cell in matrix.
  500.    * From Gary Perlman's unixstat.
  501.    */
  502.   private static double chiCell(double freq, double expect, boolean yates){
  503.     double  diff = freq - expect;
  504.        
  505.     if (yates) {
  506.       diff = Math.abs (diff) - 0.5;
  507.       if (diff < 0.0) { // over-correction
  508. diff = 0.0;
  509.       }
  510.     }
  511.     if (Math.abs(expect) < 10e-10) {
  512.       return (0.0);
  513.     } else {
  514.       return (diff * diff / expect);
  515.     }
  516.   }
  517.   /**
  518.    * Main method for testing this class.
  519.    */
  520.   public static void main(String[] ops) {
  521.     double[] firstRow = {10, 5, 20};
  522.     double[] secondRow = {2, 10, 6};
  523.     double[] thirdRow = {5, 10, 10};
  524.     double[][] matrix = new double[3][0];
  525.     matrix[0] = firstRow; matrix[1] = secondRow; matrix[2] = thirdRow;
  526.     for (int i = 0; i < matrix.length; i++) {
  527.       for (int j = 0; j < matrix[i].length; j++) {
  528. System.out.print(matrix[i][j] + " ");
  529.       }
  530.       System.out.println();
  531.     }
  532.     System.out.println("Chi-squared probability: " +
  533.        ContingencyTables.chiSquared(matrix, false));
  534.     System.out.println("Chi-squared value: " +
  535.        ContingencyTables.chiVal(matrix, false));
  536.     System.out.println("Cochran's criterion fullfilled: " +
  537.        ContingencyTables.cochransCriterion(matrix));
  538.     System.out.println("Cramer's V: " +
  539.        ContingencyTables.CramersV(matrix));
  540.     System.out.println("Entropy of first row: " +
  541.        ContingencyTables.entropy(firstRow));
  542.     System.out.println("Entropy conditioned on columns: " +
  543.        ContingencyTables.entropyConditionedOnColumns(matrix));
  544.     System.out.println("Entropy conditioned on rows: " +
  545.        ContingencyTables.entropyConditionedOnRows(matrix));
  546.     System.out.println("Entropy conditioned on rows (with Laplace): " +
  547.        ContingencyTables.entropyConditionedOnRows(matrix, matrix, 3));
  548.     System.out.println("Entropy of rows: " +
  549.        ContingencyTables.entropyOverRows(matrix));
  550.     System.out.println("Entropy of columns: " +
  551.        ContingencyTables.entropyOverColumns(matrix));
  552.     System.out.println("Gain ratio: " +
  553.        ContingencyTables.gainRatio(matrix));
  554.     System.out.println("Negative log2 of multiple hypergeometric probability: " +
  555.        ContingencyTables.log2MultipleHypergeometric(matrix));
  556.     System.out.println("Symmetrical uncertainty: " +
  557.        ContingencyTables.symmetricalUncertainty(matrix));
  558.     System.out.println("Tau value: " +
  559.        ContingencyTables.tauVal(matrix));
  560.     double[][] newMatrix = new double[3][3];
  561.     newMatrix[0][0] = 1; newMatrix[0][1] = 0; newMatrix[0][2] = 1;
  562.     newMatrix[1][0] = 0; newMatrix[1][1] = 0; newMatrix[1][2] = 0;
  563.     newMatrix[2][0] = 1; newMatrix[2][1] = 0; newMatrix[2][2] = 1;
  564.     System.out.println("Matrix with empty row and column: ");
  565.     for (int i = 0; i < newMatrix.length; i++) {
  566.       for (int j = 0; j < newMatrix[i].length; j++) {
  567. System.out.print(newMatrix[i][j] + " ");
  568.       }
  569.       System.out.println();
  570.     }
  571.     System.out.println("Reduced matrix: ");
  572.     newMatrix = ContingencyTables.reduceMatrix(newMatrix);
  573.     for (int i = 0; i < newMatrix.length; i++) {
  574.       for (int j = 0; j < newMatrix[i].length; j++) {
  575. System.out.print(newMatrix[i][j] + " ");
  576.       }
  577.       System.out.println();
  578.     }
  579.   }
  580. }