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
– Troy Rockwood, Commented Dec 13, 2013 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
– Simd, Commented Dec 13, 2013 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
– Troy Rockwood, Commented Dec 13, 2013 at 20:21

Hooked · Accepted Answer · 2013-12-13 21:38:50Z

3

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

answered Dec 13, 2013 at 21:38

Hooked

88.9k46 gold badges197 silver badges271 bronze badges

Sign up to request clarification or add additional context in comments.

1 Comment

Simd Over a year ago

Thank you. I also edited the formula slightly to make it more sensible (see the new factor of 4).

Collectives™ on Stack Overflow

Plot function with large binomial coefficients

1 Answer 1

1 Comment

Your Answer

Linked

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

1 Comment

Your Answer

Sign up or log in

Post as a guest

Linked

Related