Problem 5.

In class you used a Monte Carlo simulation to determine the probability of obtaining 60 or more heads in a trial of 100 flips of a fair coin. Let’s use the variable N for thenumber of trials used in a simulation. Modify your code so that it stores the data of the number of heads obtained in each “experiment” in an array.

(a) Let N= 100 and a plot histogram of the data using the command

plt.hist(data, 101, [30,70]])

where data is the name of your array, 101 is the number of bins, and 30 and 70 are the lowest highest values of the number of heads that you want to consider.

(b) Repeat for N = 100,000.

In [1]:
import numpy as np
from scipy import stats

import matplotlib as mpl
import matplotlib.pyplot as plt
In [2]:
# Following is an Ipython magic command that puts figures in notebook.
%matplotlib notebook
        
# M.L. modifications of matplotlib defaults
# Changes can also be put in matplotlibrc file, 
# or effected using mpl.rcParams[]
mpl.style.use('classic') 
plt.rc('figure', figsize = (6, 4.5)) # Reduces overall size of figures
plt.rc('axes', labelsize=16, titlesize=14)
plt.rc('figure', autolayout = True) # Adjusts supblot parameters for new size

(a)

In [3]:
n_flips = 100
flips=stats.randint.rvs(0,2,size=n_flips)
In [4]:
flips
Out[4]:
array([1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 0,
       1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 1, 0, 0, 0, 1, 1, 0,
       0, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0,
       0, 1, 1, 1, 1, 1, 0, 1, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1,
       0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0])
In [5]:
n_heads = np.sum(flips)
n_heads
Out[5]:
46

Now, let's do a bunch of experiments, each with 100 flips, so that we can figure what the probabilities are. E.g., prob(60 or more heads out of 100)

In [6]:
N = 100
results = np.zeros(N)
In [7]:
for i in range(N):
    flips=stats.randint.rvs(0,2,size=n_flips)
    nheads = np.sum(flips)
    results[i] = nheads
In [8]:
plt.figure()
nbins = 101
low = 30
high = 70
plt.xlabel("Number of heads out of 100")
plt.ylabel("occurences")
plt.title("Histogram; equal sized bins")
out = plt.hist(results, nbins, [low,high])
In [9]:
# What is the prob of getting 60 or more heads?
Morethan60 = 0
for i in range(N):
    if results[i] >= 60:
        Morethan60 = Morethan60 + 1
print("The prob of getting 60 or more heads out of 100 flips is ",Morethan60/N)

The prob of getting 60 or more heads out of 100 flips is  0.04

(b)

In [10]:
N = 100000
results = np.zeros(N)
In [11]:
for i in range(N):
    flips=stats.randint.rvs(0,2,size=n_flips)
    nheads = np.sum(flips)
    results[i] = nheads
In [12]:
plt.figure()
nbins = 101
low = 30
high = 70
plt.xlabel("Number of heads out of 100")
plt.ylabel("occurences")
plt.title("Histogram; equal sized bins")
out = plt.hist(results, nbins, [low,high])
In [13]:
# What is the prob of getting 60 or more heads?
Morethan60 = 0
for i in range(N):
    if results[i] >= 60:
        Morethan60 = Morethan60 + 1
print("The prob of getting 60 or more heads out of 100 flips is ",Morethan60/N)

The prob of getting 60 or more heads out of 100 flips is  0.02836