Random Numbers and Fancy Indexing

Module 3: NumPy

Random Numbers and Fancy Indexing

The Monty Hall Problem

Suppose you are on a game show, and you are given the choice of three doors: Behind one door is a car; behind the others, goats. You pick a door, say No. 1, and the host, who knows what is behind the doors, opens another door, say No. 3, which has a goat. He then says to you, “Do you want to pick door No. 2?” Is it to your advantage to switch your choice?

To answer this question, we can simulate the game a big number of times. In the first run of simulations we keep the door we chose, in the second run we switch doors. Then we compare our winning percentage.

Setup

For this lesson we will only need NumPy.

import numpy as np

ModuleNotFoundError: No module named 'numpy'

Preparing the Game

Initializing the Matrix

We can use numpy arrays to represent the doors on each game, in the form of a matrix $\boldsymbol{x}$ with shape $(N, 3)$ where $N$ is the number of games we play. Values $x=0$ represent a goat, while $x=1$ represents the car.

# Let us start with a low number of games
n_games = 10

# Initialize the matrix, filled with zeros (only goats)
x = np.zeros((n_games, 3))

print(x)

NameError: name 'np' is not defined

What is the shape of x?

print("Shape of the x:", x.shape)

NameError: name 'x' is not defined

What is the type of variable x?

print("Type of x:", type(x))

NameError: name 'x' is not defined

What about the type of its elements?

print("Types inside the x:", x.dtype)

NameError: name 'x' is not defined

By default, NumPy initializes arrays of zeros with a floating-point data type. For this game we only need 0 and 1, so we can request integers directly.

x = np.zeros((n_games, 3), dtype=int)

print("Types inside x:", x.dtype)

print(x)

NameError: name 'np' is not defined

Lets review some basic concepts:

Show the first row of the array.
Show the first column.
Show the first element.

Knowing how to do this will help in the next steps.

print("First row of x (doors of the first game):   ", x[0, :])

NameError: name 'x' is not defined

print("First column of x (first door of each game):", x[:, 0])

NameError: name 'x' is not defined

print("First element of x (first door of first game): ", x[0, 0])

NameError: name 'x' is not defined

As of now, all doors have a goat behind (a 0).

NumPy arrays are mutable, so we can assign new values to their elements.

# Add a car to the first door of the first game
x[0, 0] = 1

print(x)

NameError: name 'x' is not defined

# Add a car to the first door of every game
x[:, 0] = 1
# This is the same as saying
# For each row `:` in `x`
# Change the first value (index `0`) to `1`

print(x)

NameError: name 'x' is not defined

Now there is a car for each game, but all cars are always at the first door.

Advanced Indexing

If we want the car to be placed randomly, we must use advanced indexing.

Numpy arrays can be indexed using lists, as long as all lists have the same length.

list_0 = [a0, b0, c0]
list_1 = [a1, b1, c1]

# Being x a 2D array
x[list_0, list_1]

# Is the same as indexing the following elements
x[a0, a1], x[b0, b1], x[c0, c1]

# Sample matrix
y = np.array([
    [1, 2, 3],
    [4, 5, 6],
    [7, 8, 9]
])

# Take elements (0, 0), (1, 1) and (2, 2)
idx = [0, 1, 2]
idy = [0, 1, 2]

print(y[idx, idy])

NameError: name 'np' is not defined

If the index arrays do not have compatible shapes, indexing fails.

# Sample matrix
y = np.array([
    [1, 2, 3],
    [4, 5, 6],
    [7, 8, 9]
])

# Take elements (0, 0), (1, 1) and (2, 2)
idx = [0, 1]
idy = [0, 1, 2]

print(y[idx, idy])

Pay attention to the structure of advanced indexing:

# Good
x[[a0, b0, c0], [a1, b1, c1]]

# Bad
x[[a0, a1], [b0, b1], [c0, c1]]

# Sample matrix
y = np.array([
    [1, 2, 3],
    [4, 5, 6],
    [7, 8, 9]
])

# Wrong indexing
print(y[(0, 0), (1, 1), (2, 2)])

# Sample matrix
y = np.array([
    [1, 2, 3],
    [4, 5, 6],
    [7, 8, 9]
])

# Take elements (0, 0), (1, 1) and (2, 2)
idx = [0, 1, 2]
idy = [0, 1, 2]

# List cannot be negative
print(y[-idx, -idy])

# Sample matrix
y = np.array([
    [1, 2, 3],
    [4, 5, 6],
    [7, 8, 9]
])

# Take elements (0, 0), (1, 1) and (2, 2)
idx = [-1, -2, -3]
idy = [0, 1, 2]

# List cannot be negative, but indices inside can!
print(y[idx, idy])

NameError: name 'np' is not defined

Random Integer Generator

The advanced indexing will allow us to choose one index per row, and fill it with “1” (a car).

We just need to generate these indices randomly. For this we can use np.random module.

# Generate a random integer in the range [0, 2)
n = np.random.randint(0, 3)

print(n)

# You could import the subpackage, but it is not necessary

import numpy.random as rn

print(rn.randint(0, 10))

ModuleNotFoundError: No module named 'numpy'

# Generate 10 random integers in the range [0, 2)
n = np.random.randint(0, 3, size=10)

print(n)

# Setting a random seed ensures reproducibility of results
np.random.seed(14)  # You can use any integer here

# Generate an array of 10 random integers between 0 and 99
n = np.random.randint(0, 100, size=10)
# With the seed set, this will produce the same array every time you run it

print(n)

NameError: name 'np' is not defined

# Different seeds produce different outcomes
np.random.seed(20)
n = np.random.randint(0, 100, size=10)
print(n)

NameError: name 'np' is not defined

Knowing this information, lets create a function that generates any number $N$ of Monty Hall games, represented by a matrix $N x 3$:

def monty_hall_init(n_games: int, verbose: bool = False) -> np.ndarray:
    """Create one randomly positioned car for each simulated game."""
    games = np.zeros((n_games, 3), dtype=int)
    game_indices = np.arange(n_games)
    car_indices = np.random.randint(0, 3, size=n_games)

    if verbose:
        print("Game indices:", game_indices)
        print("Car indices: ", car_indices)

    games[game_indices, car_indices] = 1
    return games

x = monty_hall_init(1000, verbose=True)
print(x)

NameError: name 'np' is not defined

Review time! To make sure the array we generated is right, we want to check:

The minimum and maximum values, which should be 0 and 1 respectively.
The total number of cars. Should be the same as the number of games (10).
The average number of cars in total. Should be 33%.
The average number of cars per door.
The number of cars per game. Should be 1 in all cases.

print("Range of values: ", x.min(), "to", x.max())

NameError: name 'x' is not defined

print("Total number of cars:", x.sum())

NameError: name 'x' is not defined

print("Average of cars in total:", x.mean())

NameError: name 'x' is not defined

print("Average of cars per door:  ", x.mean(axis=0))

NameError: name 'x' is not defined

print("Number of cars per attempt:", x.sum(axis=1))

NameError: name 'x' is not defined

Playing the Game

Time to play the Monty Hall game!

First Attempt: Keep the Initial Pick

In this first run, we (the guest) will pick one door and stick to it. There is no need to simulate what the host reveals.

# Reinitialize the matrix of doors, this time with 10,000 games
game = monty_hall_init(10000)

print(game.mean(axis=0))

NameError: name 'monty_hall_init' is not defined

# Reinitialize the matrix of doors, this time with 10,000 games
game = monty_hall_init(10000)

# Get the number of games
n_games, n_doors = game.shape

# Initialize the number of times we find the car, as 0
sum_car = 0

# Loop through every game
for id_game in range(n_games):
  # Make a random choice between the doors (0, 1, 2)
  id_guest = np.random.randint(0, 3)
  # If we found the car, `sum_car` will increase by 1
  sum_car += game[id_game, id_guest]

print("Win probability:", sum_car / n_games)

NameError: name 'monty_hall_init' is not defined

Second Attempt: Switching Doors

To handle the host intervention, we will use sets:

There are three doors: {1, 2, 3}
Car is in door A.
Guest chooses door B, in 33% of cases B=A.
Host must reveal a door in the subset {1, 2, 3} - {A, B}. This subset may have one or two elements. We name the revealed door C.
Guest swaps to the only door available in subset {1, 2, 3} - {B, C}.

game = [0, 0, 1]
my_pick = input("Pick a door:  ")  # Im picking door 1 (index 0)
my_pick = int(my_pick)
# How can presenter pick door 3 (index 2)?

doors = {0, 1, 2}

# Find the index of the value "1" in a list
idx_car = game.index(1)
doors = doors - {idx_car, my_pick}

print(doors)

game = [0, 0, 0]
idx_car = np.random.randint(0, 3)
game[idx_car] = 1
print(game)

my_pick = input("Pick a door:  ")  # Im picking door 1 (index 0)
my_pick = int(my_pick)
# How can presenter pick door 3 (index 2)?

doors = {0, 1, 2}

# Find the index of the value "1" in a list
idx_car = game.index(1)
doors = doors - {idx_car, my_pick}

print(doors)

# Lets start with just one game
game = monty_hall_init(1)

print(game)

NameError: name 'monty_hall_init' is not defined

We guess one door at random. How can we make the host reveal one of the other two doors that have a goat?

# Make a random choice between the doors (0, 1, 2)
id_guest = np.random.randint(0, 3)
print(f"Our choice: {id_guest}")

NameError: name 'np' is not defined

np.argmax(x) returns the indices of the maximum values along an axis.

# Where is the car? Check our array with car indexes
id_car = np.argmax(game[0])

print(f"The car is here: {id_car}")

NameError: name 'np' is not defined

The host must reveal a door that is neither what we picked, nor the door with the car. This is a subset operation!

# The host can open the door that:
# - Has not been choseen by the guest
# - Does not have the car
ls_id_host = list(set([0, 1, 2]) - set([id_guest, id_car]))

print(f"The doors that the host can show are: {ls_id_host}")

NameError: name 'id_guest' is not defined

Actually, this set operation can be performed in NumPy using setdiff1d:

# The host can open the door that:
# - Has not been choseen by the guest
# - Does not have the car
ls_id_host = np.setdiff1d([0, 1, 2], [id_guest, id_car])

print(f"The doors that the host can show are: {ls_id_host}")

NameError: name 'np' is not defined

In order to choose between the list of options, we can use np.random.choice.

id_host = np.random.choice(ls_id_host)

print(f"The host reveals door: {id_host}")

NameError: name 'np' is not defined

Finally, we have to choose the door that was picked neither by us nor by the host.

# After the host reveals the door, we (the guest) change our choice
ls_id_swap = np.setdiff1d([0, 1, 2], [id_guest, id_host])

print(f"Our options to swap are: {ls_id_swap}")

NameError: name 'np' is not defined

id_swap = ls_id_swap[0]

print(f"We swap to: {id_swap}")

NameError: name 'ls_id_swap' is not defined

# Did we won?
is_car = game[0, id_swap]

print(is_car)

NameError: name 'game' is not defined

Now we can implement this logic inside a for loop:

# Start the game
game = monty_hall_init(10000)

# Get the number of games
n_games, n_doors = game.shape

# Initialize the number of times we find the car, as 0
sum_car = 0

# Loop through every game
for id_game in range(n_games):
  # Make a random choice between the doors (0, 1, 2)
  id_guest = np.random.randint(0, 3)

  # Where is the car? Check our array with car indexes
  id_car = np.argmax(game[id_game])

  # The host can open the door that:
  # - Has not been choseen by the guest
  # - Does not have the car
  ls_id_host = np.setdiff1d([0, 1, 2], [id_guest, id_car])
  id_host = np.random.choice(ls_id_host)

  # After the host reveals the door, we (the guest) change our choice
  ls_id_swap = np.setdiff1d([0, 1, 2], [id_guest, id_host])
  id_swap = ls_id_swap[0]

  # If we found the car, `sum_car` will increase by 1
  sum_car += game[id_game, id_swap]

print("Win probability:", sum_car / n_games)

NameError: name 'monty_hall_init' is not defined

ls = [2]

print(ls[0])

Alternative: To find which door should the host open, instead of set difference, we can apply the following logic.

(A) Assuming the guest chooses wrong the first time, we have two doors that the host cannot open, id_guest and id_host. We can see that the following operation will always get us that door:

id_host = (3 - id_guest - id_car)

Cases:

When guest picks 0 and car is in 1: (3 - 0 - 1) = 2 → 2 % 3 = 2
When guest picks 0 and car is in 2: (3 - 0 - 2) = 1 → 1 % 3 = 1
When guest picks 1 and car is in 2: (3 - 1 - 2) = 0 → 0 % 3 = 0

Why add % 3? Is our way to ensure the ID stays in the interval (0, 2). In the cases above it is not required but will be necessary for the following scenario.

(B) What if the guest first choice was right?

When guest picks 0: (3 - 0 - 0) = 3 → 3 % 3 = 0
When guest picks 1: (3 - 1 - 1) = 1 → 1 % 3 = 1
When the guest picks 2: (3 - 2 - 2) = -1 → -1 % 3 = 2

Then we can adjust by randomly adding 1 or 2 so that the host picks one of the other two doors.

if id_host == id_car:
  id_host = id_host + np.random.randint(1, 3)
  id_host = id_host % 3

# Start the game
game = monty_hall_init(10000)

# Get the number of games
n_games, n_doors = game.shape

# Initialize the number of times we find the car, as 0
sum_car = 0

# Loop through every game
for id_game in range(n_games):
  # Make a random choice between the doors (0, 1, 2)
  id_guest = np.random.randint(0, 3)

  # Where is the car? Check our array with car indexes
  id_car = np.argmax(game[id_game])

  # Choose the ID for the host
  id_host = (3 - id_guest - id_car) % 3
  if id_host == id_car:
    id_host += np.random.randint(1, 3)
    id_host %= 3

  # Apply the same logic to choose the door we swap to
  id_swap = (3 - id_guest - id_host) % 3

  # If we found the car, `sum_car` will increase by 1
  sum_car += game[id_game, id_swap]

print("Win probability:", sum_car / n_games)

NameError: name 'monty_hall_init' is not defined

A Look into the Future

The code we have above works fine, but it takes some time to compute one million simulations. We can speed up this process by being clever with NumPy.

By the end of module 3, we will see how to write faster NumPy code! Our key goal: remove for loops and replace it with NumPy broadcasting.

Summary and Reflection

Fancy indexing uses arrays of indices to select several positions at once.
Paired row and column index arrays select coordinate pairs, not a rectangular block.
Random integer upper bounds are excluded: randint(0, 3) can produce 0, 1, or 2.
A random seed makes a simulation reproducible.
Repeated simulation lets us estimate probabilities such as the chance of winning the Monty Hall game.

Explain why ten simulated games are not enough to estimate the switching probability reliably. What changes as the number of simulations increases?

Homework and Review

Work through the Session 16 homework exercises available here.

Homework

--- title: "Random Numbers and Fancy Indexing" subtitle: "Module 3: NumPy" format: html --- # Random Numbers and Fancy Indexing ## The Monty Hall Problem Suppose you are on a game show, and you are given the choice of three doors: Behind one door is a car; behind the others, goats. You pick a door, say No. 1, and the host, who knows what is behind the doors, opens another door, say No. 3, which has a goat. He then says to you, "Do you want to pick door No. 2?" Is it to your advantage to switch your choice? <img src="https://upload.wikimedia.org/wikipedia/commons/thumb/3/3f/Monty_open_door.svg/640px-Monty_open_door.svg.png" width="600"> To answer this question, we can simulate the game a big number of times. In the first run of simulations we keep the door we chose, in the second run we switch doors. Then we compare our winning percentage. --- ## Setup For this lesson we will only need NumPy. ```{python} import numpy as np ``` --- ## Preparing the Game ### Initializing the Matrix We can use `numpy` arrays to represent the doors on each game, in the form of a matrix $\boldsymbol{x}$ with shape $(N, 3)$ where $N$ is the number of games we play. Values $x=0$ represent a goat, while $x=1$ represents the car. ```{python} # Let us start with a low number of games n_games = 10 # Initialize the matrix, filled with zeros (only goats) x = np.zeros((n_games, 3)) print(x) ``` What is the shape of `x`? ```{python} print("Shape of the x:", x.shape) ``` What is the type of variable `x`? ```{python} print("Type of x:", type(x)) ``` What about the type of its elements? ```{python} print("Types inside the x:", x.dtype) ``` By default, NumPy initializes arrays of zeros with a floating-point data type. For this game we only need `0` and `1`, so we can request integers directly. ```{python} x = np.zeros((n_games, 3), dtype=int) print("Types inside x:", x.dtype) print(x) ``` Lets review some basic concepts: - Show the first row of the array. - Show the first column. - Show the first element. Knowing how to do this will help in the next steps. ```{python} print("First row of x (doors of the first game): ", x[0, :]) ``` ```{python} print("First column of x (first door of each game):", x[:, 0]) ``` ```{python} print("First element of x (first door of first game): ", x[0, 0]) ``` As of now, all doors have a goat behind (a `0`). NumPy arrays are mutable, so we can assign new values to their elements. ```{python} # Add a car to the first door of the first game x[0, 0] = 1 print(x) ``` ```{python} # Add a car to the first door of every game x[:, 0] = 1 # This is the same as saying # For each row `:` in `x` # Change the first value (index `0`) to `1` print(x) ``` Now there is a car for each game, but all cars are always at the first door. <img src="https://upload.wikimedia.org/wikipedia/commons/thumb/7/71/Monty-LeftCar.svg/1280px-Monty-LeftCar.svg.png" width="400"> ### Advanced Indexing If we want the car to be placed randomly, we must use advanced indexing. Numpy arrays can be indexed using lists, as long as all lists have the same length. ```python list_0 = [a0, b0, c0] list_1 = [a1, b1, c1] # Being x a 2D array x[list_0, list_1] # Is the same as indexing the following elements x[a0, a1], x[b0, b1], x[c0, c1] ``` ```{python} # Sample matrix y = np.array([ [1, 2, 3], [4, 5, 6], [7, 8, 9] ]) # Take elements (0, 0), (1, 1) and (2, 2) idx = [0, 1, 2] idy = [0, 1, 2] print(y[idx, idy]) ``` If the index arrays do not have compatible shapes, indexing fails. ```{python} #| eval: false # Sample matrix y = np.array([ [1, 2, 3], [4, 5, 6], [7, 8, 9] ]) # Take elements (0, 0), (1, 1) and (2, 2) idx = [0, 1] idy = [0, 1, 2] print(y[idx, idy]) ``` Pay attention to the structure of advanced indexing: ```python # Good x[[a0, b0, c0], [a1, b1, c1]] # Bad x[[a0, a1], [b0, b1], [c0, c1]] ``` ```{python} #| eval: false # Sample matrix y = np.array([ [1, 2, 3], [4, 5, 6], [7, 8, 9] ]) # Wrong indexing print(y[(0, 0), (1, 1), (2, 2)]) ``` ```{python} #| eval: false # Sample matrix y = np.array([ [1, 2, 3], [4, 5, 6], [7, 8, 9] ]) # Take elements (0, 0), (1, 1) and (2, 2) idx = [0, 1, 2] idy = [0, 1, 2] # List cannot be negative print(y[-idx, -idy]) ``` ```{python} # Sample matrix y = np.array([ [1, 2, 3], [4, 5, 6], [7, 8, 9] ]) # Take elements (0, 0), (1, 1) and (2, 2) idx = [-1, -2, -3] idy = [0, 1, 2] # List cannot be negative, but indices inside can! print(y[idx, idy]) ``` ### Random Integer Generator The advanced indexing will allow us to choose one index per row, and fill it with "1" (a car). We just need to generate these indices randomly. For this we can use `np.random` module. ```{python} #| eval: false # Generate a random integer in the range [0, 2) n = np.random.randint(0, 3) print(n) ``` ```{python} # You could import the subpackage, but it is not necessary import numpy.random as rn print(rn.randint(0, 10)) ``` ```{python} #| eval: false # Generate 10 random integers in the range [0, 2) n = np.random.randint(0, 3, size=10) print(n) ``` ```{python} # Setting a random seed ensures reproducibility of results np.random.seed(14) # You can use any integer here # Generate an array of 10 random integers between 0 and 99 n = np.random.randint(0, 100, size=10) # With the seed set, this will produce the same array every time you run it print(n) ``` ```{python} # Different seeds produce different outcomes np.random.seed(20) n = np.random.randint(0, 100, size=10) print(n) ``` Knowing this information, lets create a function that generates any number $N$ of Monty Hall games, represented by a matrix $N x 3$: ```{python} def monty_hall_init(n_games: int, verbose: bool = False) -> np.ndarray: """Create one randomly positioned car for each simulated game.""" games = np.zeros((n_games, 3), dtype=int) game_indices = np.arange(n_games) car_indices = np.random.randint(0, 3, size=n_games) if verbose: print("Game indices:", game_indices) print("Car indices: ", car_indices) games[game_indices, car_indices] = 1 return games x = monty_hall_init(1000, verbose=True) print(x) ``` Review time! To make sure the array we generated is right, we want to check: - The minimum and maximum values, which should be 0 and 1 respectively. - The total number of cars. Should be the same as the number of games (10). - The average number of cars in total. Should be 33%. - The average number of cars per door. - The number of cars per game. Should be 1 in all cases. ```{python} print("Range of values: ", x.min(), "to", x.max()) ``` ```{python} print("Total number of cars:", x.sum()) ``` ```{python} print("Average of cars in total:", x.mean()) ``` ```{python} print("Average of cars per door: ", x.mean(axis=0)) ``` ```{python} print("Number of cars per attempt:", x.sum(axis=1)) ``` --- ## Playing the Game Time to play the Monty Hall game! ### First Attempt: Keep the Initial Pick In this first run, we (the guest) will pick one door and stick to it. There is no need to simulate what the host reveals. ```{python} # Reinitialize the matrix of doors, this time with 10,000 games game = monty_hall_init(10000) print(game.mean(axis=0)) ``` ```{python} # Reinitialize the matrix of doors, this time with 10,000 games game = monty_hall_init(10000) # Get the number of games n_games, n_doors = game.shape # Initialize the number of times we find the car, as 0 sum_car = 0 # Loop through every game for id_game in range(n_games): # Make a random choice between the doors (0, 1, 2) id_guest = np.random.randint(0, 3) # If we found the car, `sum_car` will increase by 1 sum_car += game[id_game, id_guest] print("Win probability:", sum_car / n_games) ``` --- ### Second Attempt: Switching Doors <img src="https://upload.wikimedia.org/wikipedia/commons/thumb/1/12/Monty-RightCarSwitch.svg/1920px-Monty-RightCarSwitch.svg.png" width="400"> To handle the host intervention, we will use sets: - There are three doors: {1, 2, 3} - Car is in door A. - Guest chooses door B, in 33% of cases B=A. - Host must reveal a door in the subset {1, 2, 3} - {A, B}. This subset may have one or two elements. We name the revealed door C. - Guest swaps to the only door available in subset {1, 2, 3} - {B, C}. ```{python} #| eval: false game = [0, 0, 1] my_pick = input("Pick a door: ") # Im picking door 1 (index 0) my_pick = int(my_pick) # How can presenter pick door 3 (index 2)? doors = {0, 1, 2} # Find the index of the value "1" in a list idx_car = game.index(1) doors = doors - {idx_car, my_pick} print(doors) ``` ```{python} #| eval: false game = [0, 0, 0] idx_car = np.random.randint(0, 3) game[idx_car] = 1 print(game) my_pick = input("Pick a door: ") # Im picking door 1 (index 0) my_pick = int(my_pick) # How can presenter pick door 3 (index 2)? doors = {0, 1, 2} # Find the index of the value "1" in a list idx_car = game.index(1) doors = doors - {idx_car, my_pick} print(doors) ``` ```{python} # Lets start with just one game game = monty_hall_init(1) print(game) ``` We guess one door at random. How can we make the host reveal one of the other two doors that have a goat? ```{python} # Make a random choice between the doors (0, 1, 2) id_guest = np.random.randint(0, 3) print(f"Our choice: {id_guest}") ``` `np.argmax(x)` returns the indices of the maximum values along an axis. ```{python} # Where is the car? Check our array with car indexes id_car = np.argmax(game[0]) print(f"The car is here: {id_car}") ``` The host must reveal a door that is neither what we picked, nor the door with the car. This is a subset operation! ```{python} # The host can open the door that: # - Has not been choseen by the guest # - Does not have the car ls_id_host = list(set([0, 1, 2]) - set([id_guest, id_car])) print(f"The doors that the host can show are: {ls_id_host}") ``` Actually, this set operation can be performed in NumPy using `setdiff1d`: ```{python} # The host can open the door that: # - Has not been choseen by the guest # - Does not have the car ls_id_host = np.setdiff1d([0, 1, 2], [id_guest, id_car]) print(f"The doors that the host can show are: {ls_id_host}") ``` In order to choose between the list of options, we can use `np.random.choice`. ```{python} id_host = np.random.choice(ls_id_host) print(f"The host reveals door: {id_host}") ``` Finally, we have to choose the door that was picked neither by us nor by the host. ```{python} # After the host reveals the door, we (the guest) change our choice ls_id_swap = np.setdiff1d([0, 1, 2], [id_guest, id_host]) print(f"Our options to swap are: {ls_id_swap}") ``` ```{python} id_swap = ls_id_swap[0] print(f"We swap to: {id_swap}") ``` ```{python} # Did we won? is_car = game[0, id_swap] print(is_car) ``` Now we can implement this logic inside a for loop: ```{python} # Start the game game = monty_hall_init(10000) # Get the number of games n_games, n_doors = game.shape # Initialize the number of times we find the car, as 0 sum_car = 0 # Loop through every game for id_game in range(n_games): # Make a random choice between the doors (0, 1, 2) id_guest = np.random.randint(0, 3) # Where is the car? Check our array with car indexes id_car = np.argmax(game[id_game]) # The host can open the door that: # - Has not been choseen by the guest # - Does not have the car ls_id_host = np.setdiff1d([0, 1, 2], [id_guest, id_car]) id_host = np.random.choice(ls_id_host) # After the host reveals the door, we (the guest) change our choice ls_id_swap = np.setdiff1d([0, 1, 2], [id_guest, id_host]) id_swap = ls_id_swap[0] # If we found the car, `sum_car` will increase by 1 sum_car += game[id_game, id_swap] print("Win probability:", sum_car / n_games) ``` ```{python} ls = [2] print(ls[0]) ``` --- **Alternative:** To find which door should the host open, instead of set difference, we can apply the following logic. **(A)** Assuming the guest chooses wrong the first time, we have two doors that the host cannot open, `id_guest` and `id_host`. We can see that the following operation will always get us that door: ``` id_host = (3 - id_guest - id_car) ``` Cases: - When guest picks 0 and car is in 1: `(3 - 0 - 1) = 2` → `2 % 3 = 2` - When guest picks 0 and car is in 2: `(3 - 0 - 2) = 1` → `1 % 3 = 1` - When guest picks 1 and car is in 2: `(3 - 1 - 2) = 0` → `0 % 3 = 0` Why add `% 3`? Is our way to ensure the ID stays in the interval (0, 2). In the cases above it is not required but will be necessary for the following scenario. **(B)** What if the guest first choice was right? - When guest picks 0: `(3 - 0 - 0) = 3` → `3 % 3 = 0` - When guest picks 1: `(3 - 1 - 1) = 1` → `1 % 3 = 1` - When the guest picks 2: `(3 - 2 - 2) = -1` → `-1 % 3 = 2` Then we can adjust by randomly adding 1 or 2 so that the host picks one of the other two doors. ``` if id_host == id_car: id_host = id_host + np.random.randint(1, 3) id_host = id_host % 3 ``` ```{python} # Start the game game = monty_hall_init(10000) # Get the number of games n_games, n_doors = game.shape # Initialize the number of times we find the car, as 0 sum_car = 0 # Loop through every game for id_game in range(n_games): # Make a random choice between the doors (0, 1, 2) id_guest = np.random.randint(0, 3) # Where is the car? Check our array with car indexes id_car = np.argmax(game[id_game]) # Choose the ID for the host id_host = (3 - id_guest - id_car) % 3 if id_host == id_car: id_host += np.random.randint(1, 3) id_host %= 3 # Apply the same logic to choose the door we swap to id_swap = (3 - id_guest - id_host) % 3 # If we found the car, `sum_car` will increase by 1 sum_car += game[id_game, id_swap] print("Win probability:", sum_car / n_games) ``` --- ## A Look into the Future The code we have above works fine, but it takes some time to compute one million simulations. We can speed up this process by being clever with NumPy. <p align="center"> <img src="https://www.sciencenews.org/wp-content/uploads/2022/11/117122_no_50YA-Speed-of-light_feat-1030x580.jpg" width="400" height="auto" /> </p> By the end of module 3, we will see how to write faster NumPy code! Our key goal: remove for loops and replace it with NumPy broadcasting. --- ## Summary and Reflection - Fancy indexing uses arrays of indices to select several positions at once. - Paired row and column index arrays select coordinate pairs, not a rectangular block. - Random integer upper bounds are excluded: `randint(0, 3)` can produce `0`, `1`, or `2`. - A random seed makes a simulation reproducible. - Repeated simulation lets us estimate probabilities such as the chance of winning the Monty Hall game. Explain why ten simulated games are not enough to estimate the switching probability reliably. What changes as the number of simulations increases? --- # Homework and Review Work through the Session 16 homework exercises available here. [Homework](homework_numpy.qmd#session-16-random-numbers-and-fancy-indexing){.btn .btn-primary}