Lecture 30 – Randomness and Simulations

Data 6, Summer 2021

from datascience import *
import numpy as np
Table.interactive_plots()

from ipywidgets import interact, widgets
from IPython.display import display

Randomness

np.random.randint and seeds

Each time you run this cell, you will get a random integer between 1 and 10 inclusive.

np.random.randint(1, 11)

7

Each time you run this cell, you will get the number 9. Changing the seed (15) to something else may change the number you see. The seed is local to the cell only (no other cells in your notebook are affected by the np.random.seed(15) line).

np.random.seed(15)
np.random.randint(1, 11)

9

The size argument allows us to get back an array instead of a single number.

np.random.seed(12)
np.random.randint(25, size = 5)

array([11,  6, 17,  2,  3])

Notice that the first five numbers in the following array are the same as in the array above.

np.random.seed(12)
np.random.randint(25, size = 15)

array([11,  6, 17,  2,  3,  3, 12, 16, 22, 17, 20,  5, 13,  2, 11])

Note that when I save an array to a variable, its result doesn't change unless I call np.random.randint again.

nums = np.random.randint(10, size = 5)
nums

array([6, 0, 5, 8, 2])

# Will be the same as whatever is shown above, since we're not calling np.random.randint again
nums

array([6, 0, 5, 8, 2])

# Still the same, since we're not calling np.random.randint again
nums

array([6, 0, 5, 8, 2])

np.random.choice

coin = np.array(['heads', 'tails'])

This is equivalent to flipping a fair coin once.

np.random.choice(coin)

'tails'

This is equivalent to flipping a fair coin 10 times.

np.random.choice(coin, 10)

array(['tails', 'heads', 'heads', 'tails', 'tails', 'heads', 'tails',
       'tails', 'heads', 'heads'], dtype='<U5')

Simulations

Example: dice rolls

np.random.randint(1, 7, size = 100)

array([3, 1, 5, 4, 6, 2, 6, 1, 1, 4, 4, 2, 4, 1, 2, 2, 1, 5, 1, 5, 2, 4,
       5, 6, 4, 6, 4, 5, 6, 4, 6, 4, 5, 1, 3, 6, 4, 3, 5, 3, 3, 4, 5, 4,
       2, 5, 4, 6, 2, 5, 6, 2, 6, 4, 6, 4, 1, 4, 3, 2, 6, 5, 3, 4, 4, 3,
       5, 6, 3, 1, 1, 6, 5, 1, 4, 6, 2, 1, 2, 5, 2, 5, 6, 3, 2, 4, 4, 5,
       3, 5, 1, 1, 5, 5, 5, 4, 4, 2, 5, 3])

# Don't worry about the code, just play with the slider that appears after running.
w = widgets.FloatLogSlider(
    value=1000,
    base=10,
    min=0, # max exponent of base
    max=6, # min exponent of base
    step=0.2, # exponent step
    description='Log Slider'
)

def rolls_hist(scale):
    scale = int(scale)
    rolls = np.random.randint(1, 7, size = scale)
    fig = Table().with_columns('Number of Rolls', rolls).hist(density = False, 
                                                              bins = np.arange(0.5, 7.5),
                                                              title = f'Empirical Distribution of {scale} Dice Rolls',
                                                              show = False)
    display(fig)
    
interact(rolls_hist, scale=w);

Note: The term "empirical" means "from an experiment or simulation" (as opposed to theoretical).

Also note that the histogram is randomly generated each time you move the slider, so if you select 1000, move to some other number, and come back to 1000, the histogram will look different than it did before. (Setting the seed would prevent this from happening.)

Example: sum of dice rolls

np.random.randint(1, 7, 100)

array([5, 6, 4, 3, 6, 4, 2, 3, 2, 4, 2, 4, 5, 5, 3, 2, 6, 3, 5, 4, 3, 6,
       5, 1, 5, 5, 6, 5, 6, 1, 3, 1, 4, 2, 6, 6, 1, 2, 1, 4, 3, 1, 2, 1,
       1, 1, 3, 3, 2, 1, 1, 4, 6, 4, 4, 1, 6, 1, 2, 5, 6, 5, 4, 2, 4, 1,
       6, 1, 6, 5, 4, 6, 5, 2, 6, 6, 3, 2, 1, 4, 5, 1, 5, 3, 1, 3, 1, 5,
       4, 1, 3, 6, 1, 4, 1, 1, 2, 2, 5, 5])

This array is the result of simulating the act of rolling two die and adding their results, 100 times.

np.random.randint(1, 7, 100) + np.random.randint(1, 7, 100)

array([ 7,  8,  9,  8,  8,  7,  4,  8,  6,  4,  5, 10,  8, 10,  6,  8,  9,
        6,  3,  7,  3,  8,  6,  9, 10,  5,  5,  5,  7, 10,  7,  5,  6,  7,
       11,  6,  3,  8,  7,  9,  5,  4,  6,  9,  4, 11,  3,  5,  6,  7,  3,
        7,  9,  8,  9,  9,  7,  8,  9,  3,  6,  6,  7,  7,  3,  7,  7,  6,
       10,  9,  8,  7,  3,  7, 10, 10,  7,  3, 11,  5,  7,  7, 10,  8,  7,
        5,  6,  7,  7,  2,  7,  8, 10,  6, 10,  7,  7,  9,  4,  8])

Let's once again see what happens when we simulate this process a varying number of times.

# Don't worry about the code, just play with the slider that appears after running.
def sum_rolls_hist(scale):
    scale = int(scale)
    rolls = np.random.randint(1, 7, size = scale) + np.random.randint(1, 7, size = scale)
    fig = Table().with_columns('Number of Rolls', rolls).hist(density = False, 
                                                              bins = np.arange(1.5, 13.5),
                                                              title = f'Empirical Distribution of the Sum of Two Dice Rolls, Repeated {scale} Times',
                                                              show = False)
    display(fig)
    
interact(sum_rolls_hist, scale=w);

Quick Check 1

rolls = np.random.randint(1, 7, size = 5)
A = 2 * rolls
B = rolls + np.random.randint(1, 7, size = 5)

Repetition

coin

array(['heads', 'tails'], dtype='<U5')

flips = np.random.choice(coin, size = 10)
flips

array(['tails', 'heads', 'heads', 'tails', 'tails', 'heads', 'tails',
       'heads', 'heads', 'tails'], dtype='<U5')

# Array of 10 elements – True if flip was 'heads', False otherwise
flips == 'heads'

array([False,  True,  True, False, False,  True, False,  True,  True,
       False])

# Number of heads
np.count_nonzero(flips == 'heads')

5

Simulating Coin Flips

Idea:

1. Flip a coin 100 times. Write down the number of heads.
2. Repeat step 1 many times – say, 10,000 times.
3. Draw a histogram of the number of heads in each iteration.

num_heads_arr = np.array([])

for i in np.arange(10000):
    flips = np.random.choice(coin, size = 100)
    heads = np.count_nonzero(flips == 'heads')
    num_heads_arr = np.append(num_heads_arr, heads)

num_heads_arr

array([42., 56., 50., ..., 42., 49., 56.])

len(num_heads_arr)

10000

Table().with_columns('Number of Heads', num_heads_arr) \
       .hist(density = False, bins = np.arange(25.5, 76.5), title = 'Empirical Distribution of 100 Coin Flips')

The Birthday Problem

def one_sim():
    seen = []
    while True:
        x = np.random.randint(1, 366)
        if x in seen:
            return len(seen) + 1
        else:
            seen.append(x)

one_sim()

44

class_sizes = np.array([])
for i in range(10000):
    class_sizes = np.append(class_sizes, one_sim())

class_sizes

array([12.,  8., 54., ..., 32., 44., 48.])

Instead of drawing a histogram, we're going to do something different. Run the cell below to generate a line plot.

t = Table().with_columns('Number of Students', class_sizes).group(0)
t = t.with_columns('Probability of Shared Birthdays', 100 * t.column('count').cumsum() / 10000)

t.plot('Number of Students', 'Probability of Shared Birthdays')

It turns out that, in a class size of 23, the chances that two students share the same birthday is about 50%.

Extra: Conway's Game of Life

You're not responsible for any of the code here. Just run the following cells. Much of this code is taken from here.

Also, this code will not work in the HTML output; you need to actually open the notebook for this to work.

from time import sleep
from threading import Thread
from ipycanvas import RoughCanvas, hold_canvas

def life_step(x):
    """Game of life step"""
    nbrs_count = sum(np.roll(np.roll(x, i, 0), j, 1)
                     for i in (-1, 0, 1) for j in (-1, 0, 1)
                     if (i != 0 or j != 0))
    return (nbrs_count == 3) | (x & (nbrs_count == 2))

def draw(x, canvas, color='black'):
    with hold_canvas(canvas):
        canvas.clear()
        canvas.fill_style = '#FFF0C9'
        canvas.stroke_style = 'white'
        canvas.fill_rect(0, 0, canvas.size[0], canvas.size[1])

        canvas.fill_style = color
        canvas.stroke_style = color

        living_cells = np.where(x)
        
        rects_x = living_cells[0] * n_pixels
        rects_y = living_cells[1] * n_pixels

        canvas.fill_rects(rects_x, rects_y, n_pixels)
        canvas.stroke_rects(rects_x, rects_y, n_pixels)

This is perhaps the only relevant line of code.

We randomly initialize a 20 by 20 grid of 1s and 0s.

x = np.random.randint(2, size = (20, 20))

n_pixels = 15

canvas = RoughCanvas(width=x.shape[1]*n_pixels, height=x.shape[0]*n_pixels)
canvas.fill_style = '#FFF0C9'
canvas.stroke_style = 'white'
canvas.fill_rect(0, 0, canvas.size[0], canvas.size[1])

draw(x, canvas, '#5770B3')

canvas

class GameOfLife(Thread):
    def __init__(self, x, canvas):
        self.x = x
        self.canvas = canvas
        super(GameOfLife, self).__init__()

    def run(self):
        for _ in range(1_000):
            self.x = life_step(self.x)
            draw(self.x, self.canvas, '#5770B3')

            sleep(0.1)

GameOfLife(x, canvas).start()
canvas