我需要对矩阵 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 循环中被注释掉。
最佳答案
编辑:
看看 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%。
关于java - 费舍尔精确检验的算法或数学是什么?,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/6379058/