Statistics: Combinations In Python

Question

[Eric Normandeau]

I need to compute combinatorials (nCr) in Python but cannot find the function to do that in 'math', 'numyp' or 'stat' libraries. Something like a function of the type:

comb = calculate_combinations(n, r)

I need the number of possible combinations, not the actual combinations, so itertools.combinations does not interest me.

Finally, I want to avoid using factorials, as the numbers I'll be calculating the combinations for can get to big and the factorials are going to be monstruous.

Any clue what function from what library I should use?

Many thanks

Answer 1

[Michael Kuhn]

There is a comb function in Scipy which should do the right thing.

[Eric Normandeau]

Exactly what I was looking for. Thanks!

Answer 2

[Jarretinha]

You really don't need nothing more than math. The easiest way to calculate factorials is by using the gamma function as you can see at wikipedia. If you n,k values are small (say < 10) the best way is to hardcode them on a table (or use python bindings for GSL which already hardcoded them). Otherwise, gamma function (or its ln) is your best (and probably only) option.

I'm not sure how to do the same in python. This is how I calculate LARGE combinatorials (combinations, compositions, etc.) over protein sequences in C.

floor(exp(lgamma(N+1)-lgamma(K+1)-lgamma(N-K+1.0)));

This will give a exact value for your combination and I'm sure python has some floor() function. Otherwise, you always can use pygsl to access GSL special functions. Note that this formula do not make use of recursion and calculation of gamma function values are essentially an O(1) operation. So, this is both precise and fast.

[Eric Normandeau]

Hi @Jarretinha, my values are sometimes rather high, maybe above 1000 and I need an exact value. I would really like to compute the binomial probability directly with the combination and avoid doing the factorials myself for such large numbers. My guess is that any decent function computing the combinations will have a simplification to bypass calculating the 2 or 3 factorials needed in it. Cheers

[Jarretinha]

In this case, gammaln function is your only option. SciPy has both gamma and gammaln on special.

Answer 3

[Chris Miller]

Here's a function in Ruby that uses two tricks. First of all, it uses cancelling before calculating to remove many of the factorials that appear in both the numerator and denominator. Second, it addresses the huge number problem by doing the calculations in log2 space. I wrote this several years ago, and it seemed to solve my problem then. Should be pretty easy to translate to python.

#log base 2 definition
def log2(n)
  return Math.log(n) / Math.log(2)
end

#calculates factorial in log space
def log2factorial(n)
  sum = 0
  1.upto(n){|i|
    sum += log2(i)
  }
  return sum
end

# uses cancelling to do huge factorials
# without actually calculating them
# (a choose b)
def binomCoefficientLog2(a,b)
  top = 0
  0.upto(b-1){|i|
    top += log2(a-i)
  }                 
  bottom = log2factorial(b)
  return top - bottom
end

To get the final answer, raise two to the power of the result. Example:

> binomCoefficientLog2(10,2)
=> 5.49185309632968
> 2**binomCoefficientLog2(10,2)
=> 45.0

To verify that it works, using R:

> choose(10,2)
[1] 45

Answer 4

[Paulo Nuin]

Here is some perspectives about factorials in Python. If you follow the posts you will find that you can easily calculate combinatorial in Python with this

def binom(n, m):
    b = [0] * (n + 1)
    b[0] = 1
    for i in xrange(1, n + 1):
        b[i] = 1
        j = i - 1
        while j > 0:
            b[j] += b[j - 1]
            j -= 1
    return b[m]

and factorials.