Monte Carlo probability calculation in Python

[Rob Gaddi]

On Thu, 05 Feb 2015 08:20:41 -0800, Paul  Moore wrote:

> 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

You don't need the whole scipy stack, numpy will let you do everything 
you want.  The trick to working in numpy is to parallelize your problem; 
you don't do a thing a thousand times; you do it on a thousand-length 
array.  For example:

def dice(throws, per, sides=6):
    """Return an array throws long of rolls of (per)d(sides)."""
    all_dice = np.random.randint(1, sides+1, size=(throws, per))
    return all_dice.sum(axis=1)

-- 
Rob Gaddi, Highland Technology -- www.highlandtechnology.com
Email address domain is currently out of order.  See above to fix.