Monte Carlo estimation of the Hypergeometric Function

Question

I am trying to implement the algorithm described in the paper Statistical Test for the Comparison of Samples from Mutational Spectra (Adams & Skopek, 1986) DOI: 10.1016/0022-2836(87)90669-3:

$$p = \dfrac{\prod\limits_{i=1}^{N}(R_i!)\prod\limits_{j=1}^{M}(C_j!)}{T!\prod\limits_{i=1}^{N}\prod\limits_{j=1}^{M}(X_{ij}!)}$$

The math is easy to follow in this script although very large numbers will be generated. Would using the mpmath module be my best bet for accuracy? I notice that many of my long numbers have an excessive amount of 0s possibly due to rounding errors? I want to optimize this function to generate the p-value as accurately as possible. Speed is important but not the priority. I would prefer to implement the algorithm myself although I understand there may be hypergeometric probability functions in SciPy and PyMC. I have yet to find this specific algorithm implemented in a Python module.

1. How can I improve the precision of the very large numbers? Would working with log computations be ideal?
2. Are there any module hypergeometric functions available for this specific algorithm flavor?
3. Am I using the decimal module appropriately?

data.dat

dsum    csum
1   0
1   2
3   9
2   1
2   1
1   0
1   0
10  0
2   3
3   0
15  0
1   28
etc.

hyperg.py

from __future__ import division   # Division is now float not integer
import numpy as np
import pandas as pd
from math import factorial as fctr
from decimal import Decimal

def factorfunc(num):
    """Returns factorial of x"""
    return Decimal(fctr(num))

def main():
    """Return the Monte Carlo Hypergeometric probability"""
    frame = pd.read_csv(r'c:\dev\hyperg\data.dat', sep='\t', header=0)

    colsum_1 = sum(frame['dsum'])
    colsum_2 = sum(frame['csum'])
    full_sum = colsum_1 + colsum_2
    length = len(frame)

    frame.index = frame.index + 1
    frame2 = pd.DataFrame(frame.apply(np.sum, axis=1)).applymap(factorfunc)
    frame3 = frame.applymap(factorfunc)

    binner = 1
    for index in xrange(length):
        binner*= frame2.iloc[index]
    numerator = long(binner) * fctr(colsum_1) * fctr(colsum_2)

    binner = 1
    for index in xrange(length):
        first, second = frame3.iloc[index]
        binner *= int(first) * int(second)
    denominator = long(binner) * fctr(full_sum)
    return numerator / denominator

if __name__ == '__main__':
    HYPERGP = main()
    print HYPERGP

Answer 1

I am the Adams of Adams-Skopek. I want to point out a couple of things:

1. You get some truncation error when you sum or multiply a series of values in different orders. You need to estimate this and account for it.

2. There is a lower limit on sparsity of the data, you need to address this.

3. If you can get a copy of the Pagano and Halverson program (cited in Adams and Skopek), it is useful for comparing it's results with your results. It does the exact hypergeometric test very efficiently.

4. There is a PC version here. You might need to ask the authors for the code.

5. The computer department here can probably locate the original Fortran code for the original version of the Adams Skopek program and give you a copy (which I don't have).

Answer 2

Interpretation of data

The paper you cited says:

The algorithm for Monte Carlo estimation of the \$p\$ value of the hypergeometric test is as follows.

1. Calculate the hypergeometric probability of the \$N \times M\$ table representing the mutants observed in the experiments cross-classified by the \$M\$ treatments and the \$N\$ types and sites of mutation.

   The hypergeometric probability of the observed table, \$p\$, is given by the formula

   $$p = \dfrac{\prod\limits_{i=1}^{N}(R_i!)\prod\limits_{j=1}^{M}(C_j!)}{T!\prod\limits_{i=1}^{N}\prod\limits_{j=1}^{M}(X_{ij}!)},$$

   where the \$R_i\$ and \$C_j\$ values are the row and column marginal totals, \$T\$ is the total number of observed mutants, and the \$X_{ij}\$ values are frequencies of mutants in each cell.

I suggest writing a docstring that includes a pseudocode explanation:

'''
Monte Carlo estimation of p value of the hypergeometric test, using
formula from DOI 10.1016/0022-2836(87)90669-3:

     prod(r[i]! for i in 0:N) * prod(c[j]! for j in 0:M)
p = -----------------------------------------------------
       t! * prod(x[i][j]! for i in 0:N for j in 0:M)

where the r[i] and c[j] are the row and column marginal totals,
t is the total number of observed mutants, and x[i][j] are
frequencies of mutants in each cell.
'''

Then, use the same notation in your code!

I have a hard time seeing the correspondence between the formula and the variables in your function. One column sum makes sense (you have two). However, I don't see \$R_i\$ or \$T\$, and I don't understand where your \$X_{ij}\$ come from.

Coding practices

Instead of iterating to calculate binner, you should be able to use numpy.product(). Using higher-level functions such as product() helps readability by preventing your program from being bogged down in implementation details. You want the code to look as much like the proposed pseudocode docstring as possible.

There is a huge difference between finding factors and calculating the factorial of something. factorfunc() is a therefore a poor name. Furthermore, renaming math.factorial() to fctr() is pointlessly harmful. Hard-coding the filename names your program inflexible throwaway code. I suggest making it a command-line parameter, or at least change main() to be parametized, like main(r'c:\dev\hyperg\data.dat').