Groups | Search | Server Info | Keyboard shortcuts | Login | Register [http] [https] [nntp] [nntps]

Groups > comp.lang.python > #85257 > unrolled thread

Monte Carlo probability calculation in Python

Back to article view | Back to comp.lang.python

Contents

  Monte Carlo probability calculation in Python Paul  Moore <p.f.moore@gmail.com> - 2015-02-05 08:20 -0800
    Re: Monte Carlo probability calculation in Python Joel Goldstick <joel.goldstick@gmail.com> - 2015-02-05 11:27 -0500
      Re: Monte Carlo probability calculation in Python Paul  Moore <p.f.moore@gmail.com> - 2015-02-05 08:33 -0800
        Re: Monte Carlo probability calculation in Python Sturla Molden <sturla.molden@gmail.com> - 2015-02-07 12:30 +0000
    Re: Monte Carlo probability calculation in Python Rob Gaddi <rgaddi@technologyhighland.invalid> - 2015-02-05 16:56 +0000
      Re: Monte Carlo probability calculation in Python Paul  Moore <p.f.moore@gmail.com> - 2015-02-05 11:25 -0800
        Re: Monte Carlo probability calculation in Python Rob Gaddi <rgaddi@technologyhighland.invalid> - 2015-02-05 20:56 +0000
        Re: Monte Carlo probability calculation in Python Ian Kelly <ian.g.kelly@gmail.com> - 2015-02-05 14:00 -0700
          Re: Monte Carlo probability calculation in Python Paul  Moore <p.f.moore@gmail.com> - 2015-02-05 15:02 -0800
            Re: Monte Carlo probability calculation in Python Paul  Moore <p.f.moore@gmail.com> - 2015-02-06 14:54 -0800
              Re: Monte Carlo probability calculation in Python Steven D'Aprano <steve+comp.lang.python@pearwood.info> - 2015-02-07 10:49 +1100
                Re: Monte Carlo probability calculation in Python Paul  Moore <p.f.moore@gmail.com> - 2015-02-09 13:57 -0800
                  Re: Monte Carlo probability calculation in Python Paul  Moore <p.f.moore@gmail.com> - 2015-02-10 06:40 -0800

#85257 — Monte Carlo probability calculation in Python

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

[toc] | [next] | [standalone]

#85258

[Multipart message — attachments visible in raw view] — view raw

On Thu, Feb 5, 2015 at 11:20 AM, Paul Moore <p.f.moore@gmail.com> 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
> --
> https://mail.python.org/mailman/listinfo/python-list
>

have you googled "python monte carlo"?

-- 
Joel Goldstick
http://joelgoldstick.com

[toc] | [prev] | [next] | [standalone]

#85259

On Thursday, 5 February 2015 16:28:07 UTC, Joel Goldstick  wrote:
> have you googled "python monte carlo"?

Yes. And a number of other variations. None gave anything that seemed to relate. It's quite likely though that I'm simply not understanding how things like pymc (which came up in the searches) might help me, or how to convert my problem into a Monte Carlo integration problem (another topic that came up a lot, for example). So if there are specific links from such a search that match well to the problem as I described it above, I'd be really grateful for pointers.

Paul

[toc] | [prev] | [next] | [standalone]

#85326

Paul  Moore <p.f.moore@gmail.com> wrote:

> 
> Yes. And a number of other variations. None gave anything that seemed to
> relate. It's quite likely though that I'm simply not understanding how
> things like pymc (which came up in the searches) might help me, or how to
> convert my problem into a Monte Carlo integration problem (another topic
> that came up a lot, for example). So if there are specific links from
> such a search that match well to the problem as I described it above, I'd
> be really grateful for pointers.

PyMC and emcee are used for solving complicated integrals that often turn
up in Bayesian statistics. They have the same usecase as the R-like
application BUGS (WinBUGS). If you don't know what Markov Chain Monte
Carlo, Gibbs sampler or Metropolis-Hastings means you probably don't want
this. It is rather hardcore numerical mathematics for solving a particular
problem (multidimensional integrals that are not analytically tractable),
not something you would use for any kind of Monte Carlo simulation.

Sturla

[toc] | [prev] | [next] | [standalone]

#85261

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.

[toc] | [prev] | [next] | [standalone]

#85266

On Thursday, 5 February 2015 16:57:07 UTC, Rob Gaddi  wrote:
> 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)

Thanks, that's a help. I see the principle, but a couple of questions. With bigger problems (deal 52 cards into bridge hands a million times, for example) would memory become an issue? Also, how do you handle things that don't fit into the built-in numpy operations? (For example, Monopoly - roll 2 dice and take the sum, unless you roll a double, in which case reroll, but if you roll 3 doubles you fail - return NaN in that case).

As examples of a few more problems I'd use this for:

1. Roll 4 dice, and sum the largest 3.
2. Roll 9 dice, count how many different numbers you get.
3. Deal 4 hands of 13 cards from a full deck, determine the length of the longest suit anyone has.
4. Monopoly rolls (see above, 2 dice but reroll doubles up to 3 times)

All of these are reasonably easy to write as Python functions, but it becomes progressively harder (for a non-numpy expert) to know how to convert them into parallelizable numpy primitives.

I'll have a play with the approach you suggest though, and see how it handles some of the more complex problems my friend was showing off with :-) (his worst example was a horrendous expression to score Cribbage hands, that if I read his note correctly took about a minute to simulate 10,000,000 deals)

Paul

[toc] | [prev] | [next] | [standalone]

#85271

On Thu, 05 Feb 2015 11:25:42 -0800, Paul  Moore wrote:

> On Thursday, 5 February 2015 16:57:07 UTC, Rob Gaddi  wrote:
>> 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)
> 
> Thanks, that's a help. I see the principle, but a couple of questions.
> With bigger problems (deal 52 cards into bridge hands a million times,
> for example) would memory become an issue? Also, how do you handle
> things that don't fit into the built-in numpy operations? (For example,
> Monopoly - roll 2 dice and take the sum, unless you roll a double, in
> which case reroll, but if you roll 3 doubles you fail - return NaN in
> that case).
> 
> As examples of a few more problems I'd use this for:
> 
> 1. Roll 4 dice, and sum the largest 3.
> 2. Roll 9 dice, count how many different numbers you get.
> 3. Deal 4 hands of 13 cards from a full deck, determine the length of
> the longest suit anyone has.
> 4. Monopoly rolls (see above, 2 dice but reroll doubles up to 3 times)
> 
> All of these are reasonably easy to write as Python functions, but it
> becomes progressively harder (for a non-numpy expert) to know how to
> convert them into parallelizable numpy primitives.
> 
> I'll have a play with the approach you suggest though, and see how it
> handles some of the more complex problems my friend was showing off with
> :-) (his worst example was a horrendous expression to score Cribbage
> hands, that if I read his note correctly took about a minute to simulate
> 10,000,000 deals)
> 
> Paul

A lot of your examples can be dealt by the numpy fancy indexing tricks, 
which is really the heart of numpy.  A comparison on an array of numbers 
gives an array of bools, which can then be used to index the original 
array.  The Monopoly example, for instance.

  all_dice = np.random.randint(1, 7, size=(throws, 2))
  doubles = all_dice[all_dice[:,0] == all_dice[:,1]]

doubles is now a writeable view into the original all_dice array, so you 
could reroll them by saying something like.

  doubles[:,:] = np.random.randint(1, 7, size=doubles.shape)

Memory will eventually become a concern, but in that case you could 
always run it in parts and aggregate the parts.  I wouldn't sweat it 
until you get into 100s of millions of tries.  Even then, you could make 
some improvements by forcing the dtypes to be uint8, though there may be 
a performance penalty.

Numpy is fantastically powerful once you can wrap your head around it; I 
use it for a ton of things I used to use Matlab for.  The tutorial at  
http://wiki.scipy.org/Tentative_NumPy_Tutorial will go into lots of depth 
for you.


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

[toc] | [prev] | [next] | [standalone]

#85272

On Thu, Feb 5, 2015 at 12:25 PM, Paul  Moore <p.f.moore@gmail.com> wrote:
> On Thursday, 5 February 2015 16:57:07 UTC, Rob Gaddi  wrote:
>> 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)
>
> Thanks, that's a help. I see the principle, but a couple of questions. With bigger problems (deal 52 cards into bridge hands a million times, for example) would memory become an issue?

At the point memory becomes an issue you can partially roll it back
into a loop. For example, deal the bridge hands 10000 times in a loop
of 100.

> Also, how do you handle things that don't fit into the built-in numpy operations? (For example, Monopoly - roll 2 dice and take the sum, unless you roll a double, in which case reroll, but if you roll 3 doubles you fail - return NaN in that case).

Building on Rob's example:

def monopoly(throws, per=2, rerolls=3, sides=6):
    all_dice = np.random.randint(1, sides+1, size=(throws, rerolls, per))
    doubles = all_dice[...,0] == all_dice[...,1]
    three_doubles = doubles[:,0] & doubles[:,1] & doubles[:,2]
    return all_dice.sum(axis=2), doubles, three_doubles

This returns a (throws x rerolls) array of the sum of each roll, a
(throws x rerolls) array of booleans indicating whether the roll was a
double or not, and a throws-long array of booleans indicating whether
three doubles were rolled.

[toc] | [prev] | [next] | [standalone]

#85274

On Thursday, 5 February 2015 21:01:21 UTC, Ian  wrote:
> Building on Rob's example:
> 
> def monopoly(throws, per=2, rerolls=3, sides=6):
>     all_dice = np.random.randint(1, sides+1, size=(throws, rerolls, per))
>     doubles = all_dice[...,0] == all_dice[...,1]
>     three_doubles = doubles[:,0] & doubles[:,1] & doubles[:,2]
>     return all_dice.sum(axis=2), doubles, three_doubles
> 
> This returns a (throws x rerolls) array of the sum of each roll, a
> (throws x rerolls) array of booleans indicating whether the roll was a
> double or not, and a throws-long array of booleans indicating whether
> three doubles were rolled.

Nice! Thanks both of you. I start to see how to think about these things now. I'll play around with some of the other examples I've got, and see how far I get.

P

[toc] | [prev] | [next] | [standalone]

#85307

On Thursday, 5 February 2015 23:02:42 UTC, Paul  Moore  wrote:
> Nice! Thanks both of you. I start to see how to think about these
> things now. I'll play around with some of the other examples I've got,
> and see how far I get.

Just a quick status update, in case you're interested. With relatively little work (considering I started not knowing much about numpy) I managed to put together solutions for a couple of my friend's problems which ran in basically the same time as his custom C code. Very impressive!

Thanks again for the help.
Paul

[toc] | [prev] | [next] | [standalone]

#85312

Paul  Moore wrote:

> On Thursday, 5 February 2015 23:02:42 UTC, Paul  Moore  wrote:
>> Nice! Thanks both of you. I start to see how to think about these
>> things now. I'll play around with some of the other examples I've got,
>> and see how far I get.
> 
> Just a quick status update, in case you're interested. With relatively
> little work (considering I started not knowing much about numpy) I managed
> to put together solutions for a couple of my friend's problems which ran
> in basically the same time as his custom C code. Very impressive!


Very nice! Care to share the code?




-- 
Steven

[toc] | [prev] | [next] | [standalone]

#85407

On Friday, 6 February 2015 23:49:51 UTC, Steven D'Aprano  wrote:
> > Just a quick status update, in case you're interested. With relatively
> > little work (considering I started not knowing much about numpy) I managed
> > to put together solutions for a couple of my friend's problems which ran
> > in basically the same time as his custom C code. Very impressive!
> 
> Very nice! Care to share the code?

Will do. I'll have to do some tidying up (at the moment, it's in an IPython notebook on my PC, full of chunks of code trying various approaches) but I'll post it in a few days.

The C++ code isn't postable, unfortunately, as it's not open source. Also, it's not *quite* custom C++, but rather an expression you feed into his program that does this type of calculation. The program's built specifically for doing this type of problem, and runs an expression tree in a tight C++ loop (there's a small virtual function call overhead, but that's about it). More accurately, I guess, it's a custom C++ application for solving this class of problem. But regardless, I'd expected Python to only manage to be "good enough", not to equal the C++ application, so I'm still impressed!

But it's a nice problem as an introduction to numpy - hard enough to be worth the effort, and it *needs* the performance numpy gives. But simple enough that solving it is achievable.

Paul

[toc] | [prev] | [next] | [standalone]

#85449

On Monday, 9 February 2015 21:57:51 UTC, Paul  Moore  wrote:
> On Friday, 6 February 2015 23:49:51 UTC, Steven D'Aprano  wrote:
> > Very nice! Care to share the code?
> 
> Will do.

Here's the code I used for the Monopoly calculations.

import numpy as np

def monopoly(samples):
    # 2d6 x 3
    num = 2
    sides = 6
    rolls = 3
    dice = np.random.randint(1, sides + 1, size=(samples, rolls, num))

    # Doubles are if the two dice are the same
    doubles = dice[...,0] == dice[...,1]
    # Reroll if this (and all previous rolls) are doubles
    reroll = np.logical_and.accumulate(doubles, axis=1)
    # Keep the first entry, and any valid rerolls
    keep = np.insert(reroll, 0, True, axis=1)
    # The last column (all 3 are doubles) is "go to gaol"
    keep, gaol = np.split(keep, [-1], axis=-1)
    gaol = gaol[...,0]

    # Add up all the rolls and zero out any we don't want to keep
    totals = dice.sum(axis=-1) * keep
    # Remove any "go to gaol" cases
    totals = totals[~gaol,...]
    # Add up the distance moved
    totals = totals.sum(axis=-1)

    gaol_count = np.sum(gaol)

    return totals, gaol_count

samples = 1000000
totals, gaol_count = monopoly(samples)

print("Average distance moved =", totals.mean())
print("Percentage of times you go to gaol: {:%}".format(gaol_count/samples))
for i, count in enumerate(np.bincount(totals)):
    print("{:>2} {}".format(i, count))

The actual calculations are fairly simple, once you get your head round the tricks you need to do everything array-based. The part I find the hardest by far, is keeping track of the dimensions of the arrays involved, and what I need to do to match up axes, etc. For example, it took me a long time to work out what was going wrong with stripping out the "go to gaol" cases. The problem turned out to be that I needed the line "gaol = gaol[...,0]" to strip off a dimension - but it was really hard working out what was going on, particularly as numpy broadcasts the values you have, so you don't get errors, just weirdly wrong answers. (This isn't unique to numpy, I remember having the same issues when I used to work with J).

One thing it'd be nice to do would be to abstract out the array dimension for the number of samples, so that the "application" code could be written as if working on just one sample, and the driver code extended that to multiple samples. But I don't have the first idea how I'd do that without ending up looping in Python, which is what I need to avoid if I want performance to remain good.

Lots of scope for more things to learn here, then :-)

Thanks again to everyone that helped.
Paul

[toc] | [prev] | [standalone]

Back to top | Article view | comp.lang.python

csiph-web