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 和 cmets，请查看我的 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%。

Answer 2

这是概率方程（LaTeX 格式）：

在给定特定行和列总和的情况下获得实际矩阵的条件概率，由下式给出

[![\begin{equation}
\begin{split}
P &=\prod_{i=1}^r \prod_{j=1}^c \frac{n_{i.}!n_{.j}!}{n_{..}!n_{ij}}\\
 &=\frac{(n_{1.}!n_{2.}! \cdots n_{r.}!)(n_{.1}!n_{.2}! \cdots n_{.c}!)}{n_{..}!\prod_i \prod_j n_{ij}!}
\end{split} 
\end{equation}]

这是超几何概率函数的多元泛化。

如果您使用 100,000 次迭代，并且有较小的表（例如，最多 5x5），您将几乎接近真正精确测试的收敛。

Answer 3

我找到了我的问题的答案。今天早上与一位统计学家交谈后，他让我总结所有的值，看看结果如何。我发现值的总和（如预期的那样）高于 1。但是，我还发现我可以使用总和将 p 值缩放为 0

具有小于或等于 X^2 p 值的矩阵的条件概率值之和

除以

所有矩阵的所有条件概率值之和

我用 R 费希尔精确检验检查了我的答案