In [1]:

```
# Set up packages for lecture. Don't worry about understanding this code, but
# make sure to run it if you're following along.
import numpy as np
import babypandas as bpd
import pandas as pd
from matplotlib_inline.backend_inline import set_matplotlib_formats
import matplotlib.pyplot as plt
%reload_ext pandas_tutor
%set_pandas_tutor_options {'projectorMode': True}
set_matplotlib_formats("svg")
plt.style.use('fivethirtyeight')
np.set_printoptions(threshold=20, precision=2, suppress=True)
pd.set_option("display.max_rows", 7)
pd.set_option("display.max_columns", 8)
pd.set_option("display.precision", 2)
```

- Lab 4 is due
**tomorrow at 11:59PM**. - Homework 4 is due on
**Tuesday 10/25 at 11:59PM**. - The Midterm Project is due
**Tuesday 11/1 at 11:59PM.**Use pair programming 👯. See this post for clarifications. - 10+ more weekly office hours will be added starting next week!

The Midterm Exam is in one week, on **Friday, 10/28 during your assigned lecture**.

- It will be a 50 minute, on-paper, closed-notes exam. We will provide you with first 2 pages of the reference sheet.
- It will consist of multiple choice, fill-in-the-blank code, and short answer questions.
- Bring a pen/pencil/eraser and a photo ID. No scantron or blue book needed.
- No calculator, computers, notes, or other aids are allowed.
- Seating assignments and alternate details are coming next week.
**Today's material is on the midterm; next week's is not.**

Simulation.

- Example: What's the probability of getting 60 or more heads if we flip 100 coins?
- Example: The "Monty Hall" Problem.

- What is the probability of getting 60 or more heads if we flip 100 coins?

- While we
*could*calculate it by hand (and will learn how to in future courses), we can also approximate it using the computer:- Figure out how to do one experiment (i.e., flip 100 coins).
- Run the experiment a bunch of times.
- Find the proportion of experiments in which the number of heads was 60 or more.

- This is how we'll use
**simulation**– to approximate a probability through computation.- The techniques we will introduce in today's lecture will appear in almost every lecture for the remainder of the quarter!

- To simulate, we need a way to perform a random experiment on the computer (e.g. flipping a coin, rolling a die).

- A helpful function is
`np.random.choice(options)`

.- The input,
`options`

, is a list or array to choose from. - The output is a random element in
`options`

. By default, all elements are equally likely to be chosen.

- The input,

In [2]:

```
# Simulate a fair coin flip
np.random.choice(['Heads', 'Tails'])
```

Out[2]:

'Heads'

In [3]:

```
# Simulate a roll of a die
np.random.choice(np.arange(1, 7))
```

Out[3]:

5

`np.random.choice(options, n)`

will return an array of `n`

randomly selected elements from `options`

.

In [4]:

```
# Simulate 10 fair coin flips
np.random.choice(['Heads', 'Tails'], 10)
```

Out[4]:

array(['Tails', 'Heads', 'Tails', 'Tails', 'Heads', 'Tails', 'Heads', 'Heads', 'Heads', 'Tails'], dtype='<U5')

- By default,
`np.random.choice`

selects**with**replacement. - That is, after making a selection, that option is still available.
- e.g. if every time you draw a marble from a bag, you put it back.

- If an option can only be selected once, select
**without**replacement by specifying`replace=False`

.- e.g. if every time you draw a marble from a bag, you do not put it back.

In [5]:

```
# Choose three colleges to win free HDH swag
colleges = ['Revelle', 'John Muir', 'Thurgood Marshall',
'Earl Warren', 'Eleanor Roosevelt', 'Sixth', 'Seventh']
np.random.choice(colleges, 3, replace=False)
```

Out[5]:

array(['Earl Warren', 'Eleanor Roosevelt', 'Seventh'], dtype='<U17')

What is the probability of getting 60 or more heads if we flip 100 coins?

**Strategy:**

- Figure out how to do one experiment (i.e., flip 100 coins).
- Run the experiment a bunch of times.
- Find the proportion of experiments in which the number of heads was 60 or more.

- Use
`np.random.choice`

to flip 100 coins. - Use
`np.count_nonzero`

to count the number of heads.`np.count_nonzero(array)`

returns the number of entries in`array`

that are`True`

.

In [6]:

```
coins = np.random.choice(['Heads', 'Tails'], 100)
coins
```

Out[6]:

array(['Tails', 'Tails', 'Heads', ..., 'Heads', 'Heads', 'Tails'], dtype='<U5')

In [7]:

```
(coins == 'Heads').sum()
```

Out[7]:

49

In [8]:

```
np.count_nonzero(coins == 'Heads') # counts the number of Trues in sequence
```

Out[8]:

49

- Q: Why is it called
`count_nonzero`

? - A: In Python,
`True == 1`

and`False == 0`

, so counting the non-zero elements counts the number of`True`

s.

It's a good idea to do this, as it makes it easier to run the experiment repeatedly.

In [9]:

```
def coin_experiment():
coins = np.random.choice(['Heads', 'Tails'], 100)
return np.count_nonzero(coins == 'Heads')
```

In [10]:

```
coin_experiment()
```

Out[10]:

59

- How do we run the same code many times?
**Using a**`for`

-loop! - Each time we run the experiment, we'll need to store the results in an array.
- To do this, we'll use
`np.append`

!

- To do this, we'll use

In [11]:

```
head_counts = np.array([])
head_counts
```

Out[11]:

array([], dtype=float64)

In [12]:

```
head_counts = np.append(head_counts, 15)
head_counts
```

Out[12]:

array([15.])

In [13]:

```
head_counts = np.append(head_counts, 25)
head_counts
```

Out[13]:

array([15., 25.])

In [14]:

```
# Specify the number of repetitions
repetitions = 10000
# Create an empty array to store the results
head_counts = np.array([])
for i in np.arange(repetitions):
# For each repetition, run the experiment and add the result to head_counts
head_count = coin_experiment()
head_counts = np.append(head_counts, head_count)
```

In [15]:

```
len(head_counts)
```

Out[15]:

10000

In [16]:

```
head_counts
```

Out[16]:

array([46., 59., 51., ..., 44., 49., 45.])

In [17]:

```
# In how many experiments was the number of heads >= 60?
at_least_60 = np.count_nonzero(head_counts >= 60)
at_least_60
```

Out[17]:

289

In [18]:

```
# What is this as a proportion?
at_least_60 / repetitions
```

Out[18]:

0.0289

In [19]:

```
# Can also use np.mean()! Why?
np.mean(head_counts >= 60)
```

Out[19]:

0.0289

This is quite close to the true theoretical answer!

In [20]:

```
# The theoretical answer – don't worry about how or why this code works
import math
sum([math.comb(100, i) * (1 / 2) ** 100 for i in np.arange(60, 101)])
```

Out[20]:

0.028443966820490392

In [21]:

```
bpd.DataFrame().assign(
Number_of_Heads=head_counts
).plot(kind='hist', bins=np.arange(30, 70), density=True, ec='w', figsize=(10, 5));
plt.axvline(60, color='C1');
```

- This histogram describes the distribution of the number of heads in each experiment.
- Now we see another reason to use density histograms.
- Using density means that areas approximate
**probabilities**.

- Using density means that areas approximate

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 No. 2, and the host,

**who knows what’s 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. 1?”

**Question:**Is it to your advantage to switch your choice?

*(The question was originally posed in Parade magazine’s "Ask Marilyn" column. It is called the "Monty Hall problem" because Monty Hall was the host of the game show in question, "Let's Make a Deal.")*

You originally selected door #2. The host reveals door #3 to have a goat behind it. What should you do?

A. Might as well stick with door number #2; it has just as high a chance of winning as door #1. It doesn't matter whether you switch or not.

B. Switch to door number #1; it has a higher chance of winning than door #2.

- We'll use simulation to compute:
- The probability of winning if we switch.
- The probability of winning if we stay.
- This is just 1 - (probability of winning if we switch).

- Whichever strategy has the higher probability of winning is better!

Let's **simulate** the Monty Hall problem many times to **estimate** the probability of winning.

- Figure out how to simulate one game of Monty Hall.
- Play the game many times.
- Count the proportion of wins for each strategy (stay or switch).

When a contestant picks their door, there are three equally-likely outcomes:

- Car.
- Goat #1.
- Goat #2.

In [22]:

```
behind_picked_door = np.random.choice(['Car', 'Goat #1', 'Goat #2'])
behind_picked_door
```

Out[22]:

'Goat #2'

Suppose we can see what is behind their door (but the contestant can't).

- If it is a car, they will win if they stay.
- If it is a goat, they will win if they switch.

In [23]:

```
# Determine winning_strategy ('Stay' or 'Switch') based on what behind_picked_door is.
if behind_picked_door == 'Car':
winning_strategy = 'Stay'
else:
# A goat was behind the picked door. Monty will reveal the other goat.
# Switching wins:
winning_strategy = 'Switch'
```

Let's turn this into a function to make it easier to repeat:

In [24]:

```
def simulate_monty_hall():
behind_picked_door = np.random.choice(['Car', 'Goat #1', 'Goat #2'])
if behind_picked_door == 'Car':
winning_strategy = 'Stay'
else:
winning_strategy = 'Switch'
# print(behind_picked_door, 'was behind the door. Winning strategy:', winning_strategy)
return winning_strategy
```

In [25]:

```
simulate_monty_hall()
```

Out[25]:

'Switch'

We should save the winning strategies. To do so, let's use `np.append`

:

In [26]:

```
repetitions = 10000
winning_strategies = np.array([])
for i in np.arange(repetitions):
winning_strategy = simulate_monty_hall()
winning_strategies = np.append(winning_strategies, winning_strategy)
```

In [27]:

```
winning_strategies
```

Out[27]:

array(['Stay', 'Stay', 'Stay', ..., 'Stay', 'Switch', 'Stay'], dtype='<U32')

In [28]:

```
winning_strategies
```

Out[28]:

array(['Stay', 'Stay', 'Stay', ..., 'Stay', 'Switch', 'Stay'], dtype='<U32')

In [29]:

```
np.count_nonzero(winning_strategies == 'Switch')
```

Out[29]:

6696

In [30]:

```
np.count_nonzero(winning_strategies == 'Switch') / repetitions
```

Out[30]:

0.6696

In [31]:

```
np.count_nonzero(winning_strategies == 'Stay') / repetitions
```

Out[31]:

0.3304

- These are quite close to the true probabilities of winning per strategy ($\frac{2}{3}$ for switch, $\frac{1}{3}$ for stay).
**Conclusion:**it is better to switch.

- Looking back at our implementation, we kept track of the winning strategy in each experiment.

- However, all we really needed to keep track of was the
**number of experiments**in which the winning strategy was`'Switch'`

(or`'Stay'`

).

- Idea: Keep a
*tally*of the number of times the winning strategy was`'Switch'`

. That is, initialize`switch_count`

to 0, and add 1 to it each time the winning strategy is`'Switch'`

.

In [32]:

```
switch_count = 0
```

In [33]:

```
for i in np.arange(repetitions):
winning_strategy = simulate_monty_hall()
if winning_strategy == 'Switch':
switch_count = switch_count + 1
```

In [34]:

```
switch_count / repetitions
```

Out[34]:

0.6717

In [35]:

```
1 - switch_count / repetitions
```

Out[35]:

0.32830000000000004

No arrays needed! This strategy won't always work; it depends on the goal of the simulation.

- vos Savant asked the question in
*Parade*magazine. - She stated the correct answer:
*switch*. - She received over 10,000 letters in disagreement, including over 1,000 letters from people with Ph.D.s.

- Calculating probabilities is important, but can be hard!
- You'll learn plenty of formulas in future DSC classes.

- Simulation lets us find probabilities through computing instead of math.
- Many real-world scenarios are complicated.
- Simulation is much easier than math in many of these cases.

To estimate the probability of an event through simulation:

- Make a function that runs the experiment once.
- Run that function many, many times (usually 10000) with a
`for`

-loop, and save the results in an array with`np.append`

. - Compute the proportion of times the event occurs using
`np.count_nonzero`

.

- In the next class, we will start talking about sampling.
- Key idea: We want to learn something about a large population (e.g. all undergraduates at UCSD). However, it's far too difficult to survey everyone. If we collect a sample, what can we infer about the larger population?

**Today's material is on the midterm; next week's is not.**