Plot function with large binomial coefficients

Question

I would like to plot a function which involves binomial coefficients. The code I have is

#!/usr/bin/python
from __future__ import division
from scipy.special import binom
import matplotlib.pyplot as plt
import math

max = 500
ycoords = [sum([binom(n,w)*sum([binom(w,k)*(binom(w,k)/2**w)**(4*n/math.log(n)) for k in xrange(w+1)]) for w in xrange(1,n+1)]) for n in xrange(2,max)]

xcoords = range(2,max)

plt.plot(xcoords, ycoords)
plt.show()

Unfortunately this never terminates. If you reduce max to 40 say it works fine. Is there some way to plot this function?

I am also worried that scipy.special.binom might not be giving accurate answers as it works in floating point it seems.

Your code seems to work fine, I ran it as is and although it did take a significant amount of time, that's to be expected when you are working with 500 terms. Also, the plot worked fine. I'm unclear what the problem is. — Troy Rockwood, Dec 13 '13 at 19:53
@TroyRockwood I think that binom(n,w) is only accurate to 53 bits so the picture you get won't be right will it? — Simd, Dec 13 '13 at 19:59
That's a different question than "it never terminates" or "can you plot it", obviously the result will be limited by the accuracy of the variable precision that is used. Maybe a little indication about what the application of the question would be helpful. — Troy Rockwood, Dec 13 '13 at 20:21

score 3 · Accepted Answer · answered Dec 13 '13 at 21:38

You can get significant speedup by using numpy to compute the inner loop. First change max to N (since max is a builtin) and break up your function into smaller, more manageable chunks:

N = 500
X = np.arange(2,N)

def k_loop(w,n):
    K = np.arange(0, w+1)
    return (binom(w,K)*(binom(w,K)/2**w)**(float(n)/np.log(n))).sum()

def w_loop(n):
    v = [binom(n,w)*k_loop(w,n) for w in range(1,n+1)]
    return sum(v)

Y = [w_loop(n) for n in X]

Using N=300 as a test it takes 3.932s with the numpy code, but 81.645s using your old code. I didn't even time the N=500 case since your old code took so long!

It's worth pointing out that your function is basically exponential growth and can be approximated as such. You can see this in a semilogx plot:

enter image description here

Thank you. I also edited the formula slightly to make it more sensible (see the new factor of 4). — Simd, Dec 14 '13 at 08:35

Plot function with large binomial coefficients

1 Answers1

Linked