Monte Carlo estimation of the Hypergeometric Function I am trying to implement the algorithm described in the paper Statistical Test for the Comparison of Samples from Mutational Spectra (Adams &amp; 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.


How can I improve the precision of the very large numbers? Would working with log computations be ideal? 
Are there any module hypergeometric functions available for this specific algorithm flavor?
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
I am the Adams of Adams-Skopek.  I want to point out a couple of things:


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.
There is a lower limit on sparsity of the data, you need to address this.
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.
There is a PC version here. You might need to ask the authors for the code.
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).