probabilities, how do they work

prob.py

import math
from pylab import *

def fact(x):
    return math.gamma(x+1)

def choose(n,k):
    return fact(n) / (fact(k) * fact(n-k))

def hypergeom_mass_at_k(N, K, n, k):
    return choose(K, k) * choose(N-K, n-k) / choose(N, n)

def deal(n, k):
    values = [0] * n
    for i in xrange(0, k):
        values[i] = 1
    random.shuffle(values)
    return values

def draw_at_least_one(N, K, n):
    result = 0
    for i in xrange(1, K+1):
        try:
            result += hypergeom_mass_at_k(N, K, n, i)
        except ValueError:
            pass
    return result

def draw_at_least_one_monte_carlo(N, K, n, rounds=10000):
    count = 0
    for i in xrange(rounds):
        if sum(deal(N, K)[:n]) > 0: count += 1
    return count / float(rounds)