Monte Carlo probability calculation in Python

[Paul Moore]

I'm interested in prototyping a Monte Carlo type simulation algorithm in Python. The background is that a friend has written a similar program in C++, and I'm interested in seeing if I can achieve something comparable in a much better language :-)

The basic job of the program will be to simulate games of chance - so we'll have random inputs (die rolls, card draws, etc) and repeatedly simulate calculating a "result". Based on the results of the simulation, the idea is to estimate the probability of a given result.

So, to give a very specific example:

import random

def die(n, sides=6):
    total = sum(random.randint(1, sides) for i in range(n))
    return total

def simulate(n, test):
    "Run the simulation N times, returning the probability that TEST is true"
    successes = 0
    for i in range(n):
        if test():
            successes = successes + 1
    return successes/n

def check_3d6_gt_15():
    return die(3) > 15

if __name__ == '__main__':
    print(simulate(100000, check_3d6_gt_15))

Obviously, this is going to run ridiculously slowly as the number of simulations or the complexity of the calculation increases, but this gives the idea.

My immediate instinct is that somewhere in the scipy stack, there will be a module that does this sort of thing efficiently, but I don't really know where to look - my understanding of the maths involved is very much at the naive level above, so I'm not sure what terms I should be searching for, or how to frame a query.

Can anyone give me some pointers as to where I should go to find out more about this sort of task? Either more general theory (that would help me ask the right questions!) or specific packages or techniques I should be using in Python/numpy would be fantastic. 

Any help would be gratefully accepted - surely nobody wants to see Python beaten by a C++ program??? :-)

Thanks,
Paul