Monte Carlo probability calculation in Python
Rob Gaddi
rgaddi at technologyhighland.invalid
Thu Feb 5 11:56:19 EST 2015
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.
More information about the Python-list
mailing list