Consider the following typical probability scenario:

I defined this function to handle that scenario, I'm curious if Python has a more efficient method to handle this, or if this is the best way:

from scipy.special import binom

def grade_of_service(n, p, c):
 prob = 0
 k = c + 1
 while k <= n:
 prob += binom(n, k)*(p**k)*((1-p)**(n-k))
 k += 1
 return prob

EDIT: Here's an example solution from the author of this example: "If n = 100, p=0.1, and c=15, the probability of interest turns out to be 0.0399."
My probability does round to this result. But due to the roundoff error / numerical precision, it doesn't seem to matter if k = n or k = n+1.

Pure Python

1. You should use a for loop rather than a while loop.


2. You can use a generator comprehension to build all prob's.


3. You can use sum to get the sum.

def grade_of_service(n, p, c):
 return sum(binom(n, k)*(p**k)*((1-p)**(n-k)) for k in range(c+1, n+1))

Numpy and friends

1. Use numpy.arange rather than range.


2. Write out the equation same as above, just not in a comprehension.


3. Change sum to numpy.sum.

import numpy


def grade_of_service(n, p, c):
 k = numpy.arange(c+1, n+1)
 return numpy.sum(binom(n, k)*(p**k)*((1-p)**(n-k)))

This has a problem with numbers that exceed a certain size, as numpy numbers are finite.