Каков алгоритм ИЛИ математика для точного теста Фишера?

Question

Мне нужен точный тест Фишера для матрицы n x m. Я искал часы и нашел только один пример кода, но он написан на фортране. Я работал над Вольфрамом, и я близок к завершению, но мне не хватает самого последнего кусочка.

    /**
     * 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? Все мои попытки либо в последнем цикле for, либо в комментариях к последнему циклу for.

Вот весь код: http://pastie.org/private/f8lga9oj6f8vrxiw348q

Answer 1

Вот уравнение вероятности (в формате 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 2

edit:

Глядя на wolfram, кажется, что проблема размера nxm может быть решена с помощью:

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 и комментариев, посмотрите мой гист .

Факториалы растут очень быстро , например:

• long может хранить 20 первых факторных значений.
• doubleможет хранить 170 первых факторных значений .

Пример Wolfram:

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

приведет к:

0.023809523809523809523809523809523809523809523809523809523810

edit 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 также имеет статистические пакеты, в которых реализован точный тест Фишера.ИМХО, написание этого на Java может быть на 20% быстрее, но необходимые усилия составляют около 200%, а время разработки - 400%.

Answer 3

Я нашел ответ на свой вопрос. После разговора со статистиком сегодня утром он попросил меня подвести итог всех значений и посмотреть, что из этого получится. Я обнаружил, что сумма значений (как и ожидалось) была выше 1. Однако я также обнаружил, что могу использовать сумму для масштабирования значения p до 0.

сумма значений условной вероятности матриц с меньшими или равными X ^ 2 p-значениями

РАЗДЕЛЕНО

сумма всех значений условной вероятности всех матриц

Я проверил свой ответ с помощью точного теста Р Фишера