import numpy as np
some_numbers = np.array([3, 1, 4, 1, 5, 9, 2, 6])10 Two puzzles and more tools
10.1 Introduction
In the next chapter we will deal with some more involved problems in probability, as a preparation for statistics, where we use reasoning from probability to draw conclusions about a world like our own, where variation often appears to be more or less random.
Before we get down to the business of complex probabilistic problems in the next few chapters, let’s consider a couple of peculiar puzzles. These puzzles allow us to introduce some more of the key tools in Python for Monte Carlo resampling, and show the power of such simulation to help solve, and then reason about, problems in probability.
But, before we get to the puzzles, we will need some more elements of Python.
10.2 Selecting elements from arrays with slicing
Note 10.1: Notebook: Selecting elements by slicing
As you saw in Section 6.6, we do indexing when we put square brackets following a value that is a container, such as an array. Inside the square brackets we put another value to specify which elements we want to fetch from the container.
We will use the some_numbers array as our container for indexing:
In the indexing expression below, we have the array some_numbers, followed by the value 3 inside the square brackets, telling Python that we want the fourth value from the some_numbers container.
# Indexing with an integer inside the square brackets.
some_numbers[3]np.int64(1)In the example above, we put an integer (a whole number) inside the square brackets, to specify the position of the element we want to fetch from the container.
We can also put something called a slice inside the square brackets.
A slice specifies a range of elements to fetch from the container.
We can form a slice with an integer, followed by a colon (:), followed by another integer. The first integer specifies the start position; this is the position of the first element we want. There follows a colon. Read the colon as “up to, but not including”. Finally, we have an integer that gives the stop position. The slice, thus specified, asks Python to give us all the elements from (including) the start position, up to, but not including the stop position.
For example, here we index with a slice having start of 1 (offset 1 from the start, the position of the second element). The stop is 5, meaning we should go up to, but not include the element at position 5 (the sixth element). The result is another array, that has the elements of some_numbers from positions 1 through 4:
# Indexing with a slice (an expression including a colon).
some_numbers[1:5]In fact, we can omit the value before the colon (the start value), and Python will assume we mean 0. This indexing expression fetches the elements at position 0 through 3:
# Indexing with a slice that omits the first (start) value. Python assumes 0.
some_numbers[:4]The expression above then means “get all the elements up to (not including) position 4”, or equivalently, “get the first four elements of the array”.
We can also omit the stop value. Python assumes we mean one past the last position in the array. This is also the len of the array. By taking this stop value, the slice selects the elements starting at the start element, through to (including) the last element of the array:
# Indexing with a slice that omits the last (stop) value. Python assumes
# we mean one past the last position in the array.
some_numbers[3:]The slice in the example above therefore means “get all the elements from position 3 to the end of the array”.
Actually, we can even omit the start and the stop values, leaving just the colon. As you might expect, Python assumes 0 as the start, and one past the end as the stop, so the colon on its own means “return all the elements in the array”. As usual with slicing, we get a new array with the chosen elements.
# Just the colun, meaning "all the elements in the array".
some_numbers[:]
End of notebook: Selecting elements by slicing
selecting_by_slicing starts at Note 10.1.
10.3 The treasure fleet recovered
As promised, we have now arrived at the first of the probability puzzles.
This is a classic problem in probability:1
A Spanish treasure fleet of three ships was sunk at sea off Mexico. One ship had a chest of gold forward and another aft, another ship had a chest of gold forward and a chest of silver aft, while a third ship had a chest of silver forward and another chest of silver aft. Divers just found one of the ships and a chest of gold in it, but they don’t know whether it was from forward or aft. They are now taking bets about whether the other chest found on the same ship will contain silver or gold. What are fair odds?
These are the logical steps one may distinguish in arriving at a correct answer with deductive logic (portrayed in Figure 10.1).
Postulate three ships — Ship I with two gold chests (G-G), ship II with one gold and one silver chest (G-S), and ship III with S-S. (Choosing notation might well be considered one or more additional steps.)
Assert equal probabilities of each ship being found.
Step 2 implies equal probabilities of being found for each of the six chests.
Fact: Diver finds a chest of gold.
Step 4 implies that S-S ship III was not found; hence remove it from subsequent analysis.
Three possibilities: 6a) Diver found chest I-Ga, 6b) diver found I-Gb, 6c) diver found II-Gc.
From step 2, the cases a, b, and c in step 6 have equal probabilities.
If possibility 6a is the case, then the other chest is I-Gb; the comparable statements for cases 6b and 6c are I-Ga and II-S.
From steps 6 and 7: From equal probabilities of the three cases, and no other possible outcome, \(P(6a) = 1/3\), \(P(6b) = 1/3\), \(P(6c) = 1/3\).
So \(P(G) = P(6a) + P(6b)\) = 1/3 + 1/3 = 2/3.
See Figure 10.1.
The following simulation arrives at the correct answer.
- Write “Gold” on three pieces of paper and “Silver” on three pieces of paper. These represent the chests.
- Get three buckets each with two pieces of paper. Each bucket represents a ship, each piece of paper represents a chest in that ship. One bucket has two pieces of paper with “Gold” written on them; one has pieces of paper with “Gold” and “Silver”, and one has “Silver” and “Silver”.
- Choose a bucket at random, to represent choosing a ship at random.
- Shuffle the pieces of paper in the bucket and pick one, to represent choosing the first chest from that ship at random.
- If the piece of paper says “Silver”, the first chest we found in this ship was silver, and we stop the trial and make no further record. If “Gold”, continue.
- Get the second piece of paper from the bucket, representing the second chest on the chosen ship. Record whether this was “Silver” or “Gold” on the scoreboard.
- Repeat steps (3 - 6) many times, and calculate the proportion of “Gold”s on the scoreboard. (The answer should be about \(\frac{2}{3}\).)
Here is a notebook simulation with Python:
Note 10.2: Notebook: Ships with gold and silver
import numpy as np
rnd = np.random.default_rng()# The 3 buckets. Each bucket represents a ship. Each has two chests.
bucket1 = ['Gold', 'Gold'] # Chests in first ship.
bucket2 = ['Gold', 'Silver'] # Chests in second ship.
bucket3 = ['Silver', 'Silver'] # Chests in third ship.# For each trial, we will have one of three states:
#
# 1. When opening the first chest, it did not contain gold.
# We will reject these trials, since they do not match our
# experiment description.
# 2. Gold was found in the first and the second chest.
# 3. Gold was found in the first, but silver in the second chest.
#
# We need a placeholder value for all trials, and will make that
# "No gold in chest 1, chest 2 never opened".
second_chests = np.repeat(['No gold in chest 1, chest 2 never opened'], 10000)
for i in range(10000):
# Select a ship at random from the three ships.
ship_no = rnd.choice([1, 2, 3])
# Get the chests from this ship (represented by a bucket).
if ship_no == 1:
bucket = bucket1
if ship_no == 2:
bucket = bucket2
if ship_no == 3:
bucket = bucket3
# We shuffle the order of the chests in this ship, to simulate
# the fact that we don't know which of the two chests we have
# found first, forward or aft.
shuffled = rnd.permuted(bucket)
if shuffled[0] == 'Gold': # We found a gold chest first.
# Store whether the Second chest was silver or gold.
second_chests[i] = shuffled[1]
# End loop, go back to beginning.
# Number of times we found gold in the second chest.
n_golds = np.sum(second_chests == 'Gold')
# Number of times we found silver in the second chest.
n_silvers = np.sum(second_chests == 'Silver')
# As a ratio of golds to all second chests (where the first was gold).
print(n_golds / (n_golds + n_silvers))0.6625368731563421
End of notebook: Ships with gold and silver
gold_silver_ships starts at Note 10.2.
In the code above, we have first chosen the ship number at random, and then used a set of if ... statements to get the pair of chests corresponding to the given ship. There are simpler and more elegant ways of writing this code, but they would need some Python features that we haven’t covered yet.2
10.4 Back to Boolean arrays
The code above implements the procedure we might well use if we were simulating the problem physically. We do a trial, and we record the result. We do this on a piece of paper if we are doing a physical simulation, and in the second_chests array in code.
Finally we tally up the results. If we are doing a physical simulation, we go back over the all the trial results and counting up the “Gold” and “Silver” outcomes. In code we use the comparisons == 'Gold' and == 'Silver' to find the trials of interest, and then count them up with np.sum.
Boolean (logical) arrays are a fundamental tool in Python, and we will use them in nearly all our simulations.
Here is a remind of how those arrays work.
First, let’s slice out the first 10 values of the second_chests trial-by-trial results tally from the simulation above:
# Get values at positions 0 through 9 (up to, but not including position 10)
first_10_chests = second_chests[:10]
first_10_chestsarray(['Silver', 'No gold in chest 1, chest 2 never opened',
'No gold in chest 1, chest 2 never opened', 'Gold', 'Gold',
'No gold in chest 1, chest 2 never opened', 'Gold', 'Gold',
'No gold in chest 1, chest 2 never opened', 'Gold'], dtype='<U40')Before we started the simulation, we set second_chests to contain 10,000 strings, where each string was “No gold in chest 1, chest 2 never opened”. In the simulation, we check whether there was gold in the first chest, and, if not, we don’t change the value in second_chest, and the value remains as “No gold in chest 1, chest 2 never opened”.
Only if there was gold in the first chest, do we go on to check whether the second chest contains silver or gold. Therefore, we only set a new value in second_chests where there was gold in the first chest.
Now let’s show the effect of running a comparison on first_10_chests:
were_gold = (first_10_chests == 'Gold')
were_goldarray([False, False, False, True, True, False, True, True, False,
True])
Parentheses and Boolean comparisons
Notice the round brackets (parentheses) around (first_10_chests == 'Gold'). In this particular case, we would get the same result without the parentheses, so the paretheses are optional— although see below for an example where the they are not optional. In general, you will see we put parentheses around all expressions that generate Boolean arrays, and we recommend you do too. It is good habit to get into, to make it clear that this is an expression that generates a value.
The == 'Gold' comparison is asking a question. It is asking that question of an array, and the array contains multiple values. NumPy treats this comparison as asking the question of each element in the array. We get an answer for the question for each element. The answer for position 0 is True if the element at position 0 is equal to 'Gold' and False otherwise, and so on, for positions 1, 2 and so on. We started with 10 strings. After the comparison == 'Gold' we have 10 Boolean values, where a Boolean value can either be True or False.
Now we have an array with True for the “Gold” results and False otherwise, we can count the number of “Gold” results by using np.sum on the array. As you remember (Section 5.14) np.sum counts True as 1 and False as 0, so the sum of the Boolean array is just the number of True values in the array — the count that we need.
# The number of True values — so the number of "Gold" chests.
np.sum(were_gold)np.int64(5)10.5 Boolean arrays and another take on the ships problem
If we are doing a physical simulation, we usually want to finish up all the work for the trial during the trial, so we have one outcome from the trial. This makes it easier to tally up the results in the end.
We have no such constraint when we are using code, so it is sometimes easier to record several results from the trial, and do the final combinations and tallies at the end. We will show you what we mean with a slight variation on the two-ships code you saw above.
Note 10.3: Notebook: Another approach to ships with gold and silver
Notice that the first part of the code is identical to the first approach to this problem. There are two key differences — see the comments for an explanation.
import numpy as np
rnd = np.random.default_rng()# The 3 buckets, each representing two chests on a ship.
# As before.
bucket1 = ['Gold', 'Gold'] # Chests in first ship.
bucket2 = ['Gold', 'Silver'] # Chests in second ship.
bucket3 = ['Silver', 'Silver'] # Chests in third ship.# Here is where the difference starts. We are now going to fill in
# the result for the first chest _and_ the result for the second chest.
#
# Later we will fill in all these values, so the string we put here
# does not matter.
# Whether the first chest was Gold or Silver.
first_chests = np.repeat(['To be announced'], 10000)
# Whether the second chest was Gold or Silver.
second_chests = np.repeat(['To be announced'], 10000)
for i in range(10000):
# Select a ship at random from the three ships.
# As before.
ship_no = rnd.choice([1, 2, 3])
# Get the chests from this ship.
# As before.
if ship_no == 1:
bucket = bucket1
if ship_no == 2:
bucket = bucket2
if ship_no == 3:
bucket = bucket3
# As before.
shuffled = rnd.permuted(bucket)
# Here is the big difference - we store the result for the first and second
# chests.
first_chests[i] = shuffled[0]
second_chests[i] = shuffled[1]
# End loop, go back to beginning.
# We will do the calculation we need in the next cell. For now
# just display the first 10 values.
ten_first_chests = first_chests[:10]
print('The first 10 values of "first_chests:', ten_first_chests)The first 10 values of "first_chests: ['Gold' 'Silver' 'Silver' 'Gold' 'Gold' 'Silver' 'Gold' 'Gold' 'Silver'
'Gold']ten_second_chests = second_chests[:10]
print('The first 10 values of "second_chests', ten_second_chests)The first 10 values of "second_chests ['Silver' 'Gold' 'Silver' 'Gold' 'Gold' 'Silver' 'Gold' 'Gold' 'Silver'
'Gold']In this variant, we recorded the type of first chest for each trial (“Gold” or “Silver”), and the type of second chest of the second chest (“Gold” or “Silver”).
We would like to count the number of times there was “Gold” in the first chest and “Gold” in the second.
10.6 Combining Boolean arrays
We can do the count we need by combining the Boolean arrays with the & operator. & combines Boolean arrays with a logical and. Logical and is a rule for combining two Boolean values, where the rule is: the result is True if the first value is True and the second value if True.
Here we use the & operator to combine some Boolean values on the left and right of the operator:
True & True # Both are True, so result is TrueTrueTrue & False # At least one of the values is False, so result is FalseFalseFalse & True # At least one of the values is False, so result is FalseFalseFalse & False # At least one (in fact both) are False, result is False.False
Note 10.4:
& and and in Python
In fact Python has another operation to apply this logical and operation to values — the and operator:
print(True and True)print(True and False)print(False and True)print(False and False)You will see this and operator often in Python code, but it does not work well when combining Numpy arrays, so we will use the similar & operator, that does work on arrays.
Above you saw that the == operator (as in == 'Gold'), when applied to arrays, asks the question of every element in the array.
First make the Boolean arrays.
ten_first_gold = (ten_first_chests == 'Gold')
print("Ten first == 'Gold'", ten_first_gold)Ten first == 'Gold' [ True False False True True False True True False True]ten_second_gold = (ten_second_chests == 'Gold')
print("Ten second == 'Gold'", ten_second_gold)Ten second == 'Gold' [False True False True True False True True False True]Now let us use & to combine Boolean arrays:
ten_both = (ten_first_gold & ten_second_gold)
ten_botharray([False, False, False, True, True, False, True, True, False,
True])Notice that Python does the comparison elementwise — element by element.
You saw that when we did second_chests == 'Gold' this had the effect of asking the == 'Gold' question of each element, so there will be one answer per element in second_chests. In that case there was an array to the left of == and a single value to the right. We were comparing an array to a value.
Here we are asking the & question of ten_first_gold and ten_second_gold. Here there is an array to the left and an array to the right. We are asking the & question 10 times, but the first question we are asking is:
# First question, giving first element of result.
(ten_first_gold[0] & ten_second_gold[0])np.False_The second question is:
# Second question, giving second element of result.
(ten_first_gold[1] & ten_second_gold[1])np.False_and so on. We have ten elements on each side, and 10 answers, giving an array (ten_both) of 10 elements. Each element in ten_both is the answer to the & question for the elements at the corresponding positions in ten_first_gold and ten_second_gold.
We could also create the Boolean arrays and do the & operation all in one step, like this:
ten_both = (ten_first_chests == 'Gold') & (ten_second_chests == 'Gold')
ten_botharray([False, False, False, True, True, False, True, True, False,
True])
Parentheses, arrays and comparisons
Again you will notice the round brackets (parentheses) around (ten_first_chests == 'Gold') and (ten_second_chests == 'Gold'). Above, you saw us recommend you always use paretheses around Boolean expressions like this. The parentheses make the code easier to read — but be careful — in this case, we actually need the parentheses to make Python do what we want; see the footnote for more detail.3
Remember, we wanted the answer to the question: how many trials had “Gold” in the first chest and “Gold” in the second. We can answer that question for the first 10 trials with np.sum:
n_ten_both = np.sum(ten_both)
n_ten_bothnp.int64(5)We can answer the same question for all the trials, in the same way:
first_gold = (first_chests == 'Gold')
second_gold = (second_chests == 'Gold')
n_both_gold = np.sum(first_gold & second_gold)
n_both_goldnp.int64(3369)We could also do the same calculation all in one line:
# Notice the parentheses - we need these - see above.
n_both_gold = np.sum((first_chests == 'Gold') & (second_chests == 'Gold'))
n_both_goldnp.int64(3369)We can then count all the ships where the first chest was gold:
n_first_gold = np.sum(first_chests == 'Gold')
n_first_goldnp.int64(5085)The final calculation is the proportion of second chests that are gold, given the first chest was also gold:
p_g_given_g = n_both_gold / n_first_gold
p_g_given_gnp.float64(0.6625368731563421)Of course we won’t get exactly the same results from the two simulations, in the same way that we won’t get exactly the same results from any two runs of the same simulation, because of the random values we are using. But the logic for the two simulations are the same, and we are doing many trials (10,000), so the results will be very similar.
End of notebook: Another approach to ships with gold and silver
gold_silver_booleans starts at Note 10.3.
10.7 The Monty Hall problem
The Monty Hall Problem is a puzzle in probability that is famous for its deceptive simplicity. It has its own long Wikipedia page: https://en.wikipedia.org/wiki/Monty_Hall_problem.
Here is the problem in the form it is best known; a letter to the columnist Marilyn vos Savant, published in Parade Magazine (1990):
Suppose you’re on a game show, and you’re given the choice of three doors. Behind one door is a car, behind the others, goats. You pick a door, say #1, and the host, who knows what’s behind the doors, opens another door, say #3, which has a goat. He says to you, “Do you want to pick door #2?” Is it to your advantage to switch your choice of doors?
In fact the first person to propose (and solve) this problem was Steve Selvin, a professor of public health at the University of California, Berkeley (Selvin 1975).
Most people, including at least one of us, your humble authors, quickly come to the wrong conclusion. The most common but incorrect answer is that it will make no difference if you switch doors or stay with your original choice. The obvious intuition is that, after Monty opens his door, there are two doors that might have the car behind them, and therefore, there is a 50% chance it will be behind any one of the two. It turns out that answer is wrong; you will double your chances of winning by switching doors. Did you get the answer right?
If you got the answer wrong, you are in excellent company. As you can see from the commentary in Savant (1990), many mathematicians wrote to Parade magazine to assert that the (correct) solution was wrong. Paul Erdős was one of the most famous mathematicians of the 20th century; he could not be convinced of the correct solution until he had seen a computer simulation (Vazsonyi 1999), of the type we will do below.
To simulate a trial of this problem, we need to select a door at random to house the car, and another door at random, to be the door the contestant chooses. We number the doors 1, 2 and 3. Now we need two random choices from the options 1, 2 or 3, one for the door with the car, the other for the contestant door. To chose a door for the car, we could throw a die, and chose door 1 if the die shows 1 or 4, door 2 if the die shows 2 or 5, and door 3 for 3 or 6. Then we throw the die again to chose the contestant door.
But throwing dice is a little boring; we have to find the die, then throw it many times, and record the results. Instead we can ask the computer to chose the doors at random.
For this simulation, let us do 25 trials. We ask the computer to create two sets of 25 random numbers from 1 through 3. The first set is the door with the car behind it (“Car door”). The second set have the door that the contestant chose at random (“Our door”). We put these in a table, and make some new, empty columns to fill in later. The first new column is “Monty opens”. In due course, we will use this column to record the door that Monty Hall will open on this trial. The last two columns express the outcome. The first is “Stay wins”. This has “Yes” if we win on this trial by sticking to our original choice of door, and “No” otherwise. The last column is “Switch wins”. This has “Yes” if we win by switching doors, and “No” otherwise. See table Table 10.1).
| Car door | Our door | Monty opens | Stay wins | Switch wins | |
|---|---|---|---|---|---|
| 1 | 3 | 3 | |||
| 2 | 3 | 1 | |||
| 3 | 1 | 3 | |||
| 4 | 1 | 1 | |||
| 5 | 2 | 3 | |||
| 6 | 2 | 1 | |||
| 7 | 2 | 2 | |||
| 8 | 1 | 3 | |||
| 9 | 1 | 2 | |||
| 10 | 3 | 1 | |||
| 11 | 2 | 2 | |||
| 12 | 3 | 2 | |||
| 13 | 2 | 2 | |||
| 14 | 3 | 1 | |||
| 15 | 1 | 2 | |||
| 16 | 2 | 1 | |||
| 17 | 3 | 3 | |||
| 18 | 3 | 2 | |||
| 19 | 1 | 1 | |||
| 20 | 3 | 2 | |||
| 21 | 2 | 2 | |||
| 22 | 3 | 1 | |||
| 23 | 3 | 1 | |||
| 24 | 1 | 1 | |||
| 25 | 2 | 3 |
In the first trial in Table 10.1), the computer selected door 3 for car, and door 3 for the contestant. Now Monty must open a door, and he cannot open our door (door 3) so he has the choice of opening door 1 or door 2; he chooses randomly, and opens door 2. On this trial, we win if we stay with our original choice, and we lose if we change to the remaining door, door 1.
Now we go the second trial. The computer chose door 3 for the car, and door 1 for our choice. Monty cannot choose our door (door 1) or the door with the car behind it (door 3), so he must open door 2. Now if we stay with our original choice, we lose, but if we switch, we win.
You may want to print out table Table 10.1, and fill out the blank columns, to work through the logic.
After doing a few more trials, and some reflection, you may see that there are two different situations here: the situation when our initial guess was right, and the situation where our initial guess was wrong. When our initial guess was right, we win by staying with our original choice, but when it was wrong, we always win by switching. The chance of our initial guess being correct is 1/3 (one door out of three). So the chances of winning by staying are 1/3, and the chances of winning by switching are 2/3. But remember, you don’t need to follow this logic to get the right answer. As you will see below, the resampling simulation shows us that the Switch strategy wins.
Table 10.2 is a version of table Table 10.1 for which we have filled in the blank columns using the logic above.
| Car door | Our door | Monty opens | Stay wins | Switch wins | |
|---|---|---|---|---|---|
| 1 | 3 | 3 | 1 | Yes | No |
| 2 | 3 | 1 | 2 | No | Yes |
| 3 | 1 | 3 | 2 | No | Yes |
| 4 | 1 | 1 | 2 | Yes | No |
| 5 | 2 | 3 | 1 | No | Yes |
| 6 | 2 | 1 | 3 | No | Yes |
| 7 | 2 | 2 | 3 | Yes | No |
| 8 | 1 | 3 | 2 | No | Yes |
| 9 | 1 | 2 | 3 | No | Yes |
| 10 | 3 | 1 | 2 | No | Yes |
| 11 | 2 | 2 | 1 | Yes | No |
| 12 | 3 | 2 | 1 | No | Yes |
| 13 | 2 | 2 | 1 | Yes | No |
| 14 | 3 | 1 | 2 | No | Yes |
| 15 | 1 | 2 | 3 | No | Yes |
| 16 | 2 | 1 | 3 | No | Yes |
| 17 | 3 | 3 | 2 | Yes | No |
| 18 | 3 | 2 | 1 | No | Yes |
| 19 | 1 | 1 | 2 | Yes | No |
| 20 | 3 | 2 | 1 | No | Yes |
| 21 | 2 | 2 | 1 | Yes | No |
| 22 | 3 | 1 | 2 | No | Yes |
| 23 | 3 | 1 | 2 | No | Yes |
| 24 | 1 | 1 | 2 | Yes | No |
| 25 | 2 | 3 | 1 | No | Yes |
The proportion of times “Stay” wins in these 25 trials is 0.36. The proportion of times “Switch” wins is 0.64; the Switch strategy wins about twice as often as the Stay strategy.
10.8 Monty Hall with Python
Now you have seen what the results might look like for a physical simulation, you can exercise some of your newly-strengthened Python muscles to do the simulation with code.
Note 10.5: Notebook: The Monty Hall problem
import numpy as np
rnd = np.random.default_rng()The Monty Hall problem has a slightly complicated structure, so we will start by looking at the procedure for one trial. When we have that clear, we will put that procedure into a for loop for the simulation.
Let’s start with some variables. Let’s call the door I choose my_door.
We choose that door at random from a sequence of all possible doors. Call the doors 1, 2 and 3 from left to right.
# List of doors to chose from.
doors = [1, 2, 3]
# We choose one door at random.
my_door = rnd.choice(doors)
# Show the result
my_doornp.int64(2)We choose one of the doors to be the door with the car behind it:
# One door at random has the car behind it.
car_door = rnd.choice(doors)
# Show the result
car_doornp.int64(2)Now we need to decide which door Monty will open.
By our set up, Monty cannot open our door (my_door). By the set up, he has not opened (and cannot open) the door with the car behind it (car_door).
my_door and car_door might be the same.
So, to get Monty’s choices, we want to take all doors (doors) and remove my_door and car_door. That leaves the door or doors Monty can open.
Here are the doors Monty cannot open. Remember, a third of the time my_door and car_door will be the same, so we will include the same door twice, as doors Monty can’t open.
cant_open = [my_door, car_door]
cant_open[np.int64(2), np.int64(2)]We want to find the remaining doors from doors after removing the doors named in cant_open.
NumPy has a good function for this, called np.setdiff1d. It calculates the set difference between two sequences, such as arrays.
The set difference between two sequences is the members that are in the first sequence, but are not in the second sequence. Here are a few examples of this set difference function in NumPy.
Notice that we are using lists as the input (first and second) sequences here. We can use lists or arrays or any other type of sequence in Python. (See Section 7.3.2 for an introduction to lists).
Numpy functions like np.setdiff1d always return an array.
# Members in [1, 2, 3] that are *not* in [1]
# 1, 2, 3, removing 1, if present.
np.setdiff1d([1, 2, 3], [1])array([2, 3])# Members in [1, 2, 3] that are *not* in [2, 3]
# 1, 2, 3, removing 2 and 3, if present.
np.setdiff1d([1, 2, 3], [2, 3])array([1])# Members in [1, 2, 3] that are *not* in [2, 2]
# 1, 2, 3, removing 2 and 2 again, if present.
np.setdiff1d([1, 2, 3], [2, 2])array([1, 3])This logic allows us to choose the doors Monty can open:
montys_choices = np.setdiff1d(doors, [my_door, car_door])
montys_choicesarray([1, 3])Notice that montys_choices will only have one element left when my_door and car_door were different, but it will have two elements if my_door and car_door were the same.
Let’s play out those two cases:
my_door = 1 # For example.
car_door = 2 # For example.
# Monty can only choose door 3 now.
montys_choices = np.setdiff1d(doors, [my_door, car_door])
montys_choicesarray([3])my_door = 1 # For example.
car_door = 1 # For example.
# Monty can choose either door 2 or door 3.
montys_choices = np.setdiff1d(doors, [my_door, car_door])
montys_choicesarray([2, 3])If Monty can only choose one door, we’ll take that. Otherwise we’ll chose a door at random from the two doors available.
if len(montys_choices) == 1: # Only one door available.
montys_door = montys_choices[0] # Take the first (of 1!).
else: # Two doors to choose from:
# Choose at random.
montys_door = rnd.choice(montys_choices)
montys_doornp.int64(2)In fact, we can avoid that if len( check for the number of doors, because rnd.choice will also work on a sequence of length 1 — in that case, it always returns the single element in the sequence, like this:
# rnd.choice on sequence with single element - always returns that element.
rnd.choice([2])That means we can simplify the code above to:
# Choose single door left to choose, or door at random if two.
montys_door = rnd.choice(montys_choices)
montys_doorNow we know Monty’s door, we can identify the other door, by removing our door, and Monty’s door, from the available options:
remaining_doors = np.setdiff1d(doors, [my_door, montys_door])
# There is only one remaining door, take that.
other_door = remaining_doors[0]
other_doornp.int64(2)The logic above gives us the full procedure for one trial.
my_door = rnd.choice(doors)
car_door = rnd.choice(doors)
# Which door will Monty open?
montys_choices = np.setdiff1d(doors, [my_door, car_door])
# Choose single door left to choose, or door at random if two.
montys_door = rnd.choice(montys_choices)
# Now find the door we'll open if we switch.
remaining_doors = np.setdiff1d(doors, [my_door, montys_door])
# There is only one door left.
other_door = remaining_doors[0]
# Calculate the result of this trial.
if my_door == car_door:
stay_wins = True
if other_door == car_door:
switch_wins = TrueAll that remains is to put that trial procedure into a loop, and collect the results as we repeat the procedure many times.
# Arrays to store the results for each trial.
stay_wins = np.repeat([False], 10000)
switch_wins = np.repeat([False], 10000)
# A list of doors to chose from.
doors = [1, 2, 3]
for i in range(10000):
# You will recognize the below as the single-trial procedure above.
my_door = rnd.choice(doors)
car_door = rnd.choice(doors)
# Which door will Monty open?
montys_choices = np.setdiff1d(doors, [my_door, car_door])
# Choose single door left to choose, or door at random if two.
montys_door = rnd.choice(montys_choices)
# Now find the door we'll open if we switch.
remaining_doors = np.setdiff1d(doors, [my_door, montys_door])
# There is only one door left.
other_door = remaining_doors[0]
# Calculate the result of this trial.
if my_door == car_door:
stay_wins[i] = True
if other_door == car_door:
switch_wins[i] = True
p_for_stay = np.sum(stay_wins) / 10000
p_for_switch = np.sum(switch_wins) / 10000
print('p for stay:', p_for_stay)p for stay: 0.3326print('p for switch:', p_for_switch)p for switch: 0.6674We can also follow the same strategy as we used for the second implementation of the two-ships problem (Section 10.5).
Here, as in the second two-ships implementation, we do not calculate the trial results (stay_wins, switch_wins) in each trial. Instead, we store the doors for each trial, and then use Boolean arrays to calculate the results for all trials, at the end.
# Instead of storing the trial results, we store the doors for each trial.
my_doors = np.zeros(10000)
car_doors = np.zeros(10000)
other_doors = np.zeros(10000)
doors = [1, 2, 3]
for i in range(10000):
my_door = rnd.choice(doors)
car_door = rnd.choice(doors)
# Which door will Monty open?
montys_choices = np.setdiff1d(doors, [my_door, car_door])
# Choose single door left to choose, or door at random if two.
montys_door = rnd.choice(montys_choices)
# Now find the door we'll open if we switch.
remaining_doors = np.setdiff1d(doors, [my_door, montys_door])
# There is only one door left.
other_door = remaining_doors[0]
# Store the doors we chose.
my_doors[i] = my_door
car_doors[i] = car_door
other_doors[i] = other_door
# Now - at the end of all the trials, we use Boolean arrays to calculate the
# results.
stay_wins = my_doors == car_doors
switch_wins = other_doors == car_doors
p_for_stay = np.sum(stay_wins) / 10000
p_for_switch = np.sum(switch_wins) / 10000
print('p for stay:', p_for_stay)p for stay: 0.3374print('p for switch:', p_for_switch)p for switch: 0.662610.8.1 Insight from the Monty Hall simulation
The code simulation gives us an estimate of the right answer, but it also forces us to set out the exact mechanics of the problem. For example, by looking at the code, we see that we can calculate “stay_wins” with this code alone:
# Just choose my door and the car door for each trial.
my_doors = np.zeros(10000)
car_doors = np.zeros(10000)
doors = [1, 2, 3]
for i in range(10000):
my_doors[i] = rnd.choice(doors)
car_doors[i] = rnd.choice(doors)
# Calculate whether I won by staying.
stay_wins = my_doors == car_doors
p_for_stay = np.sum(stay_wins) / 10000
print('p for stay:', p_for_stay)p for stay: 0.3244This calculation, on its own, tells us the answer, but it also points to another insight — whatever Monty does with the doors, it doesn’t change the probability that our initial guess is right, and that must be 1 in 3 (0.333). If the probability of stay_win is 1 in 3, and we only have one other door to switch to, the probability of winning after switching must be 2 in 3 (0.666).
10.8.2 Simulation and a variant of Monty Hall
You have seen that you can avoid the silly mistakes that many of us make with probability — by asking the computer to tell you the result before you start to reason from first principles.
As an example, consider the following variant of the Monty Hall problem.
The set up to the problem has us choosing a door (my_door above), and then Monty opens one of the other two doors.
Sometimes (in fact, 2/3 of the time) there is a car behind one of Monty’s doors. We’ve obliged Monty to open the other door, and his choice is forced.
When his choice was not forced, we had Monty choose the door at random.
For example, let us say we chose door 1.
Let us say that the car is also under door 1.
Monty has the option of choosing door 2 or door 3, and he chooses randomly between them.
my_door = 1 # We chose door 1 at random.
car_door = 1 # This trial, by chance, the car door is 1.
# Monty is left with doors 2 and 3 to choose from.
montys_choices = np.setdiff1d(doors, [my_door, car_door])
# He chooses randomly.
montys_door = rnd.choice(montys_choices)
# Show the result
montys_doornp.int64(2)Now — let us say we happen to know that Monty is rather lazy, and he will always choose the left-most (lower-numbered) door of the two options.
In the previous example, Monty had the option of choosing door 2 and 3. In this new scenario, we know that he will always choose door 2 (the left-most door).
my_door = 1 # We chose door 1 at random.
car_door = 1 # This trial, by chance, the car door is 1.
# Monty is left with doors 2 and 3 to choose from.
montys_choices = np.setdiff1d(doors, [my_door, car_door])
# He chooses the left-most door, always.
montys_door = montys_choices[0]
# Show the result
montys_doornp.int64(2)It feels as if we have more information about where the car is, when we know this. Consider the situation where we have chosen door 1, and Monty opens door 3. We know that he would have preferred to open door 2, if he was allowed. We therefore know he wasn’t allowed to open door 2, and that means the car is definitely under door 2.
my_door = 1 # We chose door 1 at random.
car_door = 2 # This trial, by chance, the car door under door 2.
# Monty is left with door 3 only to choose from.
montys_choices = np.setdiff1d(doors, [my_door, car_door])
# He chooses the left-most door, always. But in this case, the left-most
# available door is 3 (he can't choose 2, it is the car_door).
# Notice the doors were in order, so the left-most door is the first door
# in the array.
montys_door = montys_choices[0]
# Show the result
montys_doornp.int64(3)To take that into account, we might try a different strategy. We will stick to our own choice if Monty has chosen the left-most of the two doors he had available to him, because he might have chosen that door because there was a car underneath the other door, or because there was a car under neither, but he preferred the left door. But, if Monty chooses the right-most of the two-doors available to him, we will switch from our own choice to the other (unopened) door, because we can be sure that the car is under the other (unopened) door.
Call this the “switch if Monty chooses right door” strategy, or “switch if right” for short.
Can you see quickly whether this will be better than the “always stay” strategy? Will it be better than the “always switch” strategy? Take a moment to think it through, and write down your answers.
If you can quickly see the answer to both questions — well done — but, are you sure you are right?
We can test by simulation.
For our test of the “switch is right” strategy, we can tell if one door is to the right of another door by comparison; higher numbers mean further to the right: 2 is right of 1, and 3 is right of 2.
# Door 3 is right of door 1.
3 > 1True# A test of the switch-if-right strategy.
# The car doors.
car_doors = np.zeros(10000)
# The door we chose using the strategy.
strategy_doors = np.zeros(10000)
doors = [1, 2, 3]
for i in range(10000):
my_door = rnd.choice(doors)
car_door = rnd.choice(doors)
# Which door will Monty open?
montys_choices = np.setdiff1d(doors, [my_door, car_door])
# Choose Monty's door from the remaining options.
# This time, he always prefers the left door.
montys_door = montys_choices[0]
# Now find the door we'll open if we switch.
remaining_doors = np.setdiff1d(doors, [my_door, montys_door])
# There is only one door remaining - but is Monty's door
# to the right of this one? Then Monty had to shift.
other_door = remaining_doors[0]
if montys_door > other_door:
# Monty's door was the right-hand door, the car is under the other one.
strategy_doors[i] = other_door
else: # We stick with the door we first thought of.
strategy_doors[i] = my_door
# Store the car door for this trial.
car_doors[i] = car_door
strategy_wins = strategy_doors == car_doors
p_for_strategy = np.sum(strategy_wins) / 10000
print('p for strategy:', p_for_strategy)p for strategy: 0.6641We find that the “switch-if-right” has around the same chance of success as the “always-switch” strategy — of about 66.6%, or 2 in 3. Were your initial answers right? Now you’ve seen the result, can you see why it should be so? It may not be obvious — the Monty Hall problem is deceptively difficult. But our case here is that the simulation first gives you an estimate of the correct answer, and then, gives you a good basis for thinking more about the problem. That is:
- simulation is useful for estimation and
- simulation is useful for reflection.
End of notebook: The Monty Hall problem
monty_hall starts at Note 10.5.
10.9 Why use simulation?
Doing these simulations has two large benefits. First, it gives us the right answer, saving us from making a mistake. Second, the process of simulation forces us to think about how the problem works. This can give us better understanding, and make it easier to reason about the solution.
We will soon see that these same advantages also apply to reasoning about statistics.
The treasure fleet problem is a restatement of a problem that Joseph Bertrand posed early in the 19th century.) Here is a variation from (Goldberg 1986, 99):
↩︎Three identical boxes each contain two coins. In one box both are pennies, in the second both are nickels, and in the third there is one penny and one nickel.
A man chooses a box at random and takes out a coin. If the coin is a penny, what is the probability that the other coin in the box is also a penny?
Although we have used multiple
ifstatements to set the chests for each ship in the code for the ship problem, it would have been more elegant to store the definitions of ships and chests as a list of lists. We haven’t covered these, but to give you a taste of what that would look like, for a single trial.↩︎import numpy as np rnd = np.random.default_rng() # We define the ships and their respective chests, using a list of lists. # Remember, list has elements, where the elements can be any type of value. # In this case, the elements in the list are themselves - lists. ship_chests = [ # Create list, with lists as elements. ['Gold', 'Gold'], # List for ship 1; element 0 of the list of lists. ['Gold', 'Silver'], # List for ship 2; element 1 of the list of lists. ['Silver', 'Silver'] # Ship 3; element 2. ] # Choose a ship at random. ship_no = rnd.choice([0, 1, 2]) # Get the corresponding chests for this ship. chests = ship_chests[ship_no] chests['Gold', 'Silver']We warned that we need parentheses around our
&expressions to get the result we want. We would add the parentheses in any case, as good practice, but here we also need the parentheses in(ten_first_chests == 'Gold') & (ten_second_chests == 'Gold'). Remember operator precedence; for example, the multiply operator*has higher precedence than the operator+, so3 + 5 * 2is equal to3 + (5 * 2)= 13. If we want to do addition before multiplication, we use parentheses to tell Python the order it should use:(3 + 5) * 2= 16.The same applies for the two operators
==and&here. In fact&has a higher precedence than==. This means that, if we write the expression without parentheses —ten_first_chests == 'Gold' & ten_second_chests == 'Gold'— because of operator precedence, Python takes this to meanten_first_chests == ('Gold' & ten_second_chests) == 'Gold'. Python does not know what to do with'Gold' & ten_second_chestsand generates an error of form'bitwise_and' not supported for the input types. The error tells you that Python does not know how to apply&('bitwise_and') to the string'Gold’ and the arrayten_second_chests.This is the same error you would get for running the code
'Gold' & ten_second_chestson its own.The point to take away is, that when you are using
&to combine Boolean arrays in Python, remember operator precedence, and, when in doubt, put parentheses around the expressions on either side of&, as here.↩︎