In my earlier post I used a simple model of the human fruit machine to calculate the expected average payout for a few variations of the game, and commented in a qualitative way on why in some cases the actual payout might vary more or less. The model actually allows us to be more specific and, for example, predict how often any given payout will occur. How well the model describes the actual *human* fruit machine is a question I will look at in a separate post.

So far I have used basic probability theory to calculate the results I wanted, and it would certainly be possible to continue that way. However, another approach is to use a computer program to simulate the fruit machine and to run it many times to see how often each outcome occurs. This doesn't give the same precision, especially for outcomes that are rare and therefore happen only a few times if at all in the simulation, but it does have some advantages:

- It is easy to change the rules of the game and see what happens without additional laborious calculations.
- It provides a different way of understanding the model, which may be more intuitive.
- Having found out how to use MathJax to display mathematical formulas in this blog last time, I now get to play with the IPython notebook. (Thanks to Daniel Rodriguez for his IPython plugin for Pelican.)

Let's start by writing a function to simulate a single round, with one fruit picked from each box. Rather than give names to the types of fruit we'll just use the numbers 0, 1, 2, 3, 4.

```
from random import Random
r = Random()
def pick_fruit(num_boxes=3, num_fruit=5):
fruit = []
for box in xrange(num_boxes):
this_fruit = r.randint(0, num_fruit - 1)
fruit.append(this_fruit)
return tuple(fruit)
```

We can test this informally by running it a few times to see if the results "look random":

```
for round in xrange(5):
print pick_fruit()
```

To simulate the flow of money over a series of games, such as all the games played at the Fun Day, we need to call this function in a loop and apply the rules regarding how much to pay out if a prize is won, and when to give the player an extra turn. Programming is easiest when each function is short and does only one thing, so let's start by writing a function to calculate the payout for a given combination of fruit:

```
def calc_prize(fruit):
""" Prize is £10 if all fruit are of type 0, or £2 for all matching of a different type. """
prize = True
first_fruit = fruit[0]
for other_fruit in fruit[1:]:
if not other_fruit == first_fruit:
prize = False
amount = 0
if prize:
if first_fruit == 0:
amount = 10
else:
amount = 2
return amount
```

Then we had better test this:

```
assert calc_prize((0, 0, 0)) == 10 # 3 bananas win £10
assert calc_prize((2, 2, 2)) == 2 # 3 of something else win £2
assert calc_prize((0, 1, 0)) == 0 # 2 the same win nothing
```

It's good I ran these tests, since the first time round I forgot to initialize `amount`

to zero in the function, resulting in an error whenever the fruit didn't all match.

We also need a function to decide whether to allow another go:

```
def another_go(fruit):
"""Returns True if any two fruit match. For simplicity only deal with the case of three boxes."""
assert len(fruit) == 3
if fruit[0] == fruit[1] and not fruit[0] == fruit[2]:
return True
if fruit[0] == fruit[2] and not fruit[0] == fruit[1]:
return True
if fruit[1] == fruit[2] and not fruit[1] == fruit[0]:
return True
return False
```

A few tests:

```
assert another_go((0, 1, 0))
assert another_go((3, 2, 2))
assert not another_go((1, 2, 3))
assert not another_go((2, 2, 2))
```

Now we are in a position to simulate a series of games:

```
def series_payout(num_games):
payout = 0
for game in xrange(num_games):
keep_going = True
while keep_going:
fruit = pick_fruit()
payout += calc_prize(fruit)
keep_going = another_go(fruit)
return payout
```

We can run this a few times to see how the payout varies. Let's take 42 as a reasonable number of games that might be played in one day, since that is the number that were actually recorded on the Fun Day.

```
num_games = 42
for i in xrange(5):
print series_payout(num_games)
```

To get a better idea of how this might work out in practice, we want to repeat this many times and plot a histogram of the results:

```
%pylab inline
num_games = 42
num_trials = 1000
income = num_games * 1 # income is £1 for each game
profits = []
for i in xrange(num_trials):
payout = series_payout(num_games)
profit = income - payout
profits.append(profit)
histo_profits = plt.hist(profits)
```

The results would seem to indicate that we are quite likely to make a profit of at least £20, but that there is some chance of making a small overall loss. Running this program again after changing the rules or the prizes is simple, and we don't have to redo demanding calculations if we want to experiment with changing the assumptions of the model, such as seeing what happens if some fruit are more likely to be picked than others. By increasing the number of games, we can also see that the small probability of making a loss vanishes almost completely for more than about 50 games.

Another example of something we can investigate using our model is what cash float we are likely to need, i.e. how much money we should have available at the beginning of the day to be reasonably sure we have taken in enough money to pay them. The float needed is simply the maximum cumulative loss at any point during the series of games, minus £1 if this loss is more than zero since we always gain an extra pound at the start of a game.

```
def series_float_needed(num_games):
cash = 0
max_loss = 0
for game in xrange(num_games):
cash += 1 # income is £1 for each game
keep_going = True
while keep_going:
fruit = pick_fruit()
cash -= calc_prize(fruit)
keep_going = another_go(fruit)
if cash < (- max_loss):
max_loss = - cash
if max_loss > 0:
return max_loss - 1
else:
return 0
```

We can simulate any number of series of games, keeping track of the float needed in each case, and plot the results as a histogram:

```
floats_needed = []
for trial in xrange(1000):
flt = series_float_needed(42)
floats_needed.append(flt)
histo_floats_needed = plt.hist(floats_needed)
```

Most of the time we would get away with not having a float, but to avoid embarrassment we should have £20 or so on hand.

Because we are almost certain to make a profit in the long term, any losses early in the game are likely to be recouped later. Thus increasing the number of games in each series beyond about 20 makes very little difference to the results. With some patience we can also increase the number of trials to look for the much rarer events where a larger loss is made, but even with 10000 trials of series of 100 games each, I did not find any cases where a float of more than £25 was needed.

Of course, if I were presenting this result in its own right rather than as an example of exploring a system using Python, I would choose my bins carefully and label my axes. But for an initial exploration, the default options at least ensure I don't overlook any high outlying values.

## Comments