Monte Carlo coin flip simulation

Question

I've been learning about Monte Carlo simulations on MIT's intro to programming class, and I'm trying to implement one that calculates the probability of flipping a coin heads side up 4 times in a row out of ten flips.

Basically, I calculate if the current flip in a 10 flip session is equal to the prior flip, and if it is, I increment a counter. Once that counter has reached 3, I exit the loop even if I haven't done all 10 coin flips since subsequent flips have no bearing on the probability. Then I increment a counter counting the number of flip sessions that successfully had 4 consecutive heads in a row. At the end, I divide the number of successful sessions by the total number of trials.

The simulation runs 10,000 trials.

def simThrows(numFlips):
    consecSuccess = 0   ## number of trials where 4 heads were flipped consecutively
    coin = 0  ## to be assigned to a number between 0 and 1
    numTrials = 10000
    for i in range(numTrials):
        consecCount = 0
        currentFlip = ""
        priorFlip = ""
        while consecCount <= 3:
            for flip in range(numFlips):
                coin = random.random()
                if coin < .5:
                    currentFlip = "Heads"
                    if currentFlip == priorFlip:
                        consecCount += 1
                        priorFlip = "Heads"
                    else:
                        consecCount = 0
                        priorFlip = "Heads"
                elif coin > .5:
                    currentFlip = "Tails"
                    if currentFlip == priorFlip:
                        consecCount = 0
                        priorFlip = "Tails"
                    else:
                        consecCount = 0
                        priorFlip = "Tails"
            break
        if consecCount >= 3:
            print("Success!")
            consecSuccess += 1
    print("Probability of getting 4 or more heads in a row: " + str(consecSuccess/numTrials))

simThrows(10)

The code seems to work (gives me consistent results each time that seem to line up to the analytic probability), but it's a bit verbose for such a simple task. Does anyone see where my code can get better?

Answer 1

When you're writing code for educational purposes (or sometimes other purposes), verbose is good because it helps you understand what's really going on. So making the code shorter or snappier or whatever is not necessarily going to make it better.

With that disclaimer out of the way: one of the most common ways to condense Python code is to use list comprehensions or generators instead of loops. A list comprehension is what you use when you're constructing a list element by element: in its simplest form, instead of this,

the_list = []
for something in something_else:
    the_list.append(func(something))

you write this:

the_list = [func(something) for something in something_else]

If you're doing something else instead of creating a list, you can have Python create an object that generates the elements on demand, rather than actually creating a list out of them. An object of that sort is called a generator and you can create one like this:

the_generator = (func(something) for something in something_else)

You can omit the parentheses when the generator is passed to another function as an argument, though.

the_sum = sum(func(something) for something in something_else)

would be equivalent to, but better than,

count = 0
for something in something_else:
    count += func(something)

There are a lot of functions in Python that take iterables (list, generators, etc.) and "condense" them into one value using some sort of operation. You can also create your own, corresponding to whatever you would be doing to the result of the loop. You can convert most loops into generator expressions this way.

So let's investigate how you could use a generator to represent the sequences of consecutive throws in each trial. You can create a generator that produces 10 random numbers easily:

(random.random() for i in xrange(10))

(this is for Python 2.x; xrange was renamed to range for Python 3). Or you can create a generator that produces 10 random values which are either 0 or 1:

(random.randint(0,1) for i in xrange(10))

That saves you from having to check each random number against 0.5. In fact, you could produce a generator that produces 10 randomly chosen words, "Heads" or "Tails", like so:

(random.choice(("Heads","Tails")) for i in xrange(10))

but it'll be easier to stick with numbers. (It's usually better to represent things with numbers or objects than with strings.)

But perhaps you're thinking, "why are you telling me to make 10 numbers when I only have to check until I find a group of four consecutive heads?" For one thing, if you're just flipping 10 coins each time, it really doesn't matter because you'll make the computer flip at most 6, and on average 3, extra coins in each trial. That doesn't take very long - it'll extend the runtime of this part of your program by 50%, but we're talking 50% of a fraction of a second. It's not worth the effort to figure out how to do it for such a small number of flips. But if each flip had, say, a billion trials, then you would definitely want to stop early. Fortunately, a generator can do this for you! Since generators produce their elements only on demand, you can stop taking elements from it once you get what you want, and not waste much of any computation. I'll address this more later.

Anyway, suppose we have our generator that produces 10 binary values 0 (tails) or 1 (heads). Is there a way to go through this and check to see whether there is a sequence of four or more heads? It turns out that just such a function is provided in itertools.groupby, which takes any iterable (list, generator, etc.) and groups consecutive identical elements. An example of its usage is

for k, g in itertools.groupby([1,0,0,1,1,1,1,0,0,0]):
    print k, list(g)

and this would print out something like

1 [1]
0 [0,0]
1 [1,1,1,1]
0 [0,0,0]

So you can check for four or more consecutive heads by just looking at the length of the group and whether the key is heads or tails.

for k, g in itertools.groupby(random.randint(0,1) for i in xrange(10)):
    if k and len(g) >= 4:
        # got a run of 4 or more consecutive heads!
        # wait, what now?

(In Python, 1 is true and 0 is false in a boolean context, so if k is equivalent to if k == 1.) OK, what shall we do with our run of 4 or more consecutive heads? Well, you're trying to find the number of trials in which this occurs. So it probably makes sense to set a success flag if this happens.

success = False
for k, g in itertools.groupby(random.randint(0,1) for i in xrange(10)):
    if k and len(g) >= 4:
        success = True
        break # this stops asking the generator for new values

But wait! This is starting to look a lot like the kind of loop that can be converted to a generator expression, isn't it? The only catch is that we're not adding anything up or constructing a list. But there is another function, any, that will go through a generator until it finds an element which matches a condition, and that's just what this for loop does. So you could write this as

success = any(k and len(g) >= 4 for k, g in
                itertools.groupby(random.randint(0,1) for i in xrange(10))

Now finally, you'll want to count how many times this happens over, say, 10000 trials. So you might write that as something like this:

successes = 0
for i in xrange(10000):
    if any(k and len(g) >= 4 for k, g in
                itertools.groupby(random.randint(0,1) for i in xrange(10)):
        successes += 1

But of course, we can also convert this to a generator, since you're just adding up numbers:

successes = sum(1 for i in xrange(10000)
      if any(k and len(g) >= 4 for k, g in
              itertools.groupby(random.randint(0,1) for i in xrange(10)))

The generator produces a 1 each time it finds a group of 4 consecutive 1s among the 10 random numbers generated.

The last thing you'd want to do is divide by the total number of trials. Well, actually, what you really want to do is calculate the average instead of the sum, and in some places you can find a function mean which is kind of like sum except that it calculates the mean instead of the total. You could use such a function if you had it. But I don't know that one is in the Python standard library, so you can just do the division:

probability = sum(1 for i in xrange(10000)
      if any(k and len(g) >= 4 for k, g in
        itertools.groupby(random.randint(0,1) for i in xrange(10))) / 10000

So the task you're trying to accomplish can actually be written in one line of Python. But it's a rather complicated line, and I wouldn't necessarily recommend actually doing this. Sometimes it's good to use a good old fashioned for loop to keep the code clear. More often, though, it's better to split your code up into modular pieces that are more useful than just what you're using them for.

Answer 2

@Keith Randall showed a short soltion, but verbosity is not necessary a bad thing. Often an important thing is to decomopose the task into separate sub-tasks, which are generic and reusable:

import random

def montecarlo(experiment, trials):
    r = 0
    for i in range(trials):
        experiment.reset()
        r += 1 if experiment.run() else 0
    return r / float(trials)


class Coinrowflipper():
    def __init__(self, trials, headsinrow):
        self.headsinrow = headsinrow
        self.trials = trials
        self.reset()
    def reset(self):
        self.count = 0
    def run(self):
        for i in range(self.trials):
            if random.random() < 0.5:
                self.count += 1
                if self.count == self.headsinrow:
                    return True
            else:
                self.count = 0
        return False

c = Coinrowflipper(10,4)
print montecarlo(c, 1000)

This way you have a generic monte-carlo experiment, and a generic heads-in a row solution.

Note: using strings ("Heads") to mark internal state is discouraged in any programming language...

Answer 3

So a few errors/simplifications I can see:

1. You're incrementing consecCount if you have either two heads or two tails.
2. setting the priorFlip value can be done regardless of what the prior flip was
3. Your while loop doesn't do anything as you have a for loop in it, then break right after it, so you'll always make all 10 tosses. A while loop isn't constantly checking the value inside of it; only when it finishes executing the code block. As you break after the for loop, the while condition will never be checked
4. There's no reason to track the previous flip, if you always reset the count on a tail
5. while the odds are small, you'll be throwing away any toss that is equal to .5

Anyway, here's my first pass at a simplified version:

def simThrows(numFlips):
    consecSuccess = 0   ## number of trials where 4 heads were flipped consecutively
    coin = 0  ## to be assigned to a number between 0 and 1
    numTrials = 10000
    for i in range(numTrials):
        consecCount = 0
        for flip in range(numFlips):
            coin = random.random()
            #since random.random() generates [0.0, 1.0), we'll make the lower bound inclusive
            if coin <= .5:
                #Treat the lower half as heads, so increment the count
                consecCount += 1 
            # No need for an elif here, it wasn't a head, so it must have been a tail
            else:
                # Tails, so we don't care about anything, just reset the count
                consecCount = 0 
            #After each flip, check if we've had 4 consecutive, otherwise keep flipping
            if consecCount >= 3:
                print("Success!")
                consecSuccess += 1
                # we had enough, so break out of the "for flip" loop
                # we'll end up in the for i in range loop . . .
                break            
    print("Probability of getting 4 or more heads in a row: " + str(consecSuccess/numTrials))

Answer 4

Based on @Karoly Horvath but different enough to be a standalone answer.

A class is very much overkill for this problem, I have achieved the same functionality with a function.

import random

def montecarlo(experiment, trials):
    positive = 0
    for _ in range(trials):
        positive += 1 if experiment() else 0
    return positive / float(trials)

def successive_heads(trials,heads_in_a_row):
    count = 0
    for _ in range(trials):
        if random.random() < 0.5:
            count += 1
            if count == heads_in_a_row:
                return True
        else:
            count = 0
    return False

print(montecarlo(lambda: successive_heads(100,10), 10000))