java - Fisher 精确检验的算法或数学是什么？

Question

我需要对矩阵 n x m 进行 Fisher 精确检验。我已经搜索了几个小时，但只找到了一个示例代码，但它是用 Fortran 语言编写的。我一直在使用 Wolfram 工作，我快要完成了，但我错过了最后一点。

    /**
     * Performs Fisher's Exact Test on a matrix m x n
     * @param matrix Any matrix m x n.
     * @return The Fisher's Exact value of the matrix
     * @throws IllegalArgumentException If the rows are not of equal length
     * @author Ryan Amos
     */
    public static double getFisherExact(int[][] matrix){
        System.out.println("Working with matrix: ");
        printMatrix(matrix);
        for (int[] array : matrix) {
            if(array.length != matrix[0].length)
                throw new IllegalArgumentException();
        }
        boolean chiSq = matrix.length != 2 || matrix[0].length != 2;
        int[] rows = new int[matrix.length];
        int[] columns = new int[matrix[0].length];
        int n;
        //compute R and C values
        for (int i = 0; i < matrix.length; i++) {
            for (int j = 0; j < matrix[i].length; j++) {
                rows[i] += matrix[i][j];
                columns[j] += matrix[i][j];
            }
            System.out.println("rows[" + i + "] = " + rows[i]);
        }

        for (int i = 0; i < columns.length; i++) {
            System.out.println("columns[" + i + "] = " + columns[i]);
        }

        //compute n
        n = 0;
        for (int i = 0; i < columns.length; i++) {
            n += columns[i];
        }

        int[][][] perms = findAllPermutations(rows, columns);
        double sum = 0;
        //int count = 0;
        double cutoff = chiSq ? getChiSquaredValue(matrix, rows, columns, n) : getConditionalProbability(matrix, rows, columns, n);
        System.out.println("P cutoff = " + cutoff + "\n");
        for (int[][] is : perms) {
            System.out.println("Matrix: ");
            printMatrix(is);
            double val = chiSq ? getChiSquaredValue(is, rows, columns, n) : getConditionalProbability(is, rows, columns, n);
            System.out.print("Value: " + val); 
            if(val <= cutoff){
                //count++;
                System.out.print(" is below " + cutoff);
//              sum += (chiSq) ? getConditionalProbability(is, rows, columns, n) : val;
//              sum += val;
                double p = getConditionalProbability(is, rows, columns, n);
                System.out.print("\np = " + p + "\nsum = " + sum + " + p = ");
                sum += p;
                System.out.print(sum);
            } else {
                System.out.println(" is above " + cutoff + "\np = " + getConditionalProbability(is, rows, columns, n));
            }
            System.out.print("\n\n");
        }
        return sum;
        //return count / (double)perms.length;
    }

所有其他方法都已经过测试和调试。问题是我不确定从哪里找到所有可能的矩阵(具有相同行和列总和的所有矩阵)。我不确定如何使用这些矩阵并将它们转换为 p 值。我读了一些关于卡方的东西，所以我找到了一个卡方算法。

所以我的问题是:根据我所拥有的(矩阵的所有排列)，我如何计算 p 值？我所有的尝试要么在最后一个 for 循环中，要么在最后一个 for 循环中被注释掉。

这是完整的代码:http://pastie.org/private/f8lga9oj6f8vrxiw348q

Answer 1

编辑:

看看 wolfram，似乎 n x m 大小的问题可以用以下方法解决:

public static BigDecimal getHypergeometricDistribution(//
        int a[][], int scale, int roundingMode//
) throws OutOfMemoryError, NullPointerException {
    ArrayList<Integer> R = new ArrayList<Integer>();
    ArrayList<Integer> C = new ArrayList<Integer>();
    ArrayList<Integer> E = new ArrayList<Integer>();
    int n = 0;

    for (int i = 0; i < a.length; i++) {
        for (int j = 0; j < a[i].length; j++) {
            if (a[i][j] < 0)
                return null;

            n += a[i][j];
            add(C, j, a[i][j]);
            add(R, i, a[i][j]);
            E.add(a[i][j]);
        }
    }
    BigDecimal term1 = //
    new BigDecimal(multiplyFactorials(C).multiply(multiplyFactorials(R)));
    BigDecimal term2 = //
    new BigDecimal(getFactorial(n).multiply(multiplyFactorials(E)));

    return term1.divide(term2, scale, roundingMode);
}

对于 getBinomialCoefficient、getFactorial 和评论，查看我的 gist .

阶乘增长非常快，例如:

• long 可以存储 20 个阶乘值。
• double can store 170 first factorial values .

Wolfram 示例案例:

    int[][] a = { { 5, 0 }, { 1, 4 } };
    System.out.println(hdMM.getHypergeometricDistribution(a, 60, 6));

会导致:

0.023809523809523809523809523809523809523809523809523809523810

编辑 2:

我的方法很快，但内存效率不高，如果输入矩阵元素的总和超过 10000，这可能是个问题。原因是阶乘的内存。

Mathematica 中几乎等价的函数，没有这个问题:

FeT1::usage = "Fisher's exact Test, 1 tailed. For more information:
    http://mathworld.wolfram.com/FishersExactTest.html";
FeT1[a_List, nr_Integer: 6] := Module[{},
   SumRow[array_] := Total[Transpose[array]]; 
   SumTotal[array_] := Total[Total[array]]; 
   SumColumn[array_] := Total[array]; 
   TF[list_] := Times @@ (list!); 
   N[(TF[SumColumn[a]]*TF[SumRow[a]])/(SumTotal[a]!* TF[Flatten[a]]), nr]
 ]; 

和示例用法:

a = {{5, 0}, {1, 4}};
FeT1[a, 59]

会屈服于

0.023809523809523809523809523809523809523809523809523809523810

Mathematica 也有可用的统计包，其中实现了 Fisher 精确检验。恕我直言，用 Java 编写这个可以快 20%，但所需的工作量大约为 200%，开发时间为 400%。