Lecture 16– Probability I¶

by Lisa Yan

In [2]:
import numpy as np
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
import seaborn as sns

# big font helper
def adjust_fontsize(size=None):
    SMALL_SIZE = 8
    MEDIUM_SIZE = 10
    BIGGER_SIZE = 12
    if size != None:
        SMALL_SIZE = MEDIUM_SIZE = BIGGER_SIZE = size

    plt.rc('font', size=SMALL_SIZE)          # controls default text sizes
    plt.rc('axes', titlesize=SMALL_SIZE)     # fontsize of the axes title
    plt.rc('axes', labelsize=MEDIUM_SIZE)    # fontsize of the x and y labels
    plt.rc('xtick', labelsize=SMALL_SIZE)    # fontsize of the tick labels
    plt.rc('ytick', labelsize=SMALL_SIZE)    # fontsize of the tick labels
    plt.rc('legend', fontsize=SMALL_SIZE)    # legend fontsize
    plt.rc('figure', titlesize=BIGGER_SIZE)  # fontsize of the figure title

plt.style.use('fivethirtyeight')
sns.set_context("talk")
sns.set_theme()
#plt.style.use('default') # revert style to default mpl
adjust_fontsize(size=20)
%matplotlib inline

Plotting helper functions.

In [34]:
# helper functions to plot and 
# compute expectation, variance, standard deviation
def plot_dist(dist_df,
                      xname="x", pname="P(X = x)", varname="X",
                      save=False):
    """
    Plot a distribution from a distribution table.
    Single-variate.
    """
    plt.bar(dist_df[xname], dist_df[pname])
    plt.ylabel(pname)
    plt.xlabel(xname)
    plt.title(f"Distribution of ${varname}$")
    plt.xticks(sorted(dist_df[xname].unique()))
    if save:
        fig = plt.gcf()
        fig.patch.set_alpha(0.0)
        plt.savefig(f"dist{varname}.png", bbox_inches = 'tight');


def simulate_samples(df, xname="x", pname="P(X = x)", size=1):
    return np.random.choice(
                df[xname], # draw from these choiecs
                size=size, # this many times
                p=df[pname]) # according to this distribution

def simulate_iid_df(dist_df, nvars, rows, varname="X"):
    """
    Make an (row x nvars) dataframe
    by calling simulate_samples for each of the nvars per row
    """
    sample_dict = {}
    for i in range(nvars):
        # generate many datapoints 
        sample_dict[f"{varname}_{i+1}"] = \
            simulate_samples(dist_df, size=rows)
    return pd.DataFrame(sample_dict)


def plot_simulated_dist(df, colname, show_stats=True, save=False, **kwargs):
    """
    Plot a simulated population.
    """
    sns.histplot(data=df, x=colname, stat='probability', discrete=True, **kwargs)
    plt.xticks(sorted(df[colname].unique())) # if there are gaps)
    if show_stats:
        display(stats_df_multi(df, [colname]))
    if save:
        fig = plt.gcf()
        fig.patch.set_alpha(0.0)
        plt.savefig(f"sim{colname}.png", bbox_inches = 'tight');

def stats_df_multi(df, colnames):
    means = df[colnames].mean(axis=0)
    variances = df[colnames].var(axis=0)
    stdevs = df[colnames].std(axis=0)
    df_stats = pd.concat([means, variances, stdevs],axis=1).T
    df_stats['index_col'] = ["E[•]", "Var(•)", "SD(•)"]
    df_stats = df_stats.set_index('index_col', drop=True).rename_axis(None)
    return df_stats

def plot_simulated_dist_multi(df, colnames, show_stats=True):
    """
    If multiple columns provided, use separate plots.
    """
    ncols = 1
    nrows = len(colnames)
    plt.figure(figsize=(6, 2*nrows+2))
    
    for i, colname in enumerate(colnames):
        subplot_int = int(100*int(nrows) + 10*int(ncols) + int(i+1))
        plt.subplot(subplot_int)
        plot_simulated_dist(df, colname, show_stats=False)
    plt.tight_layout()
    if show_stats:
        display(stats_df_multi(df, colnames))

A Random Variable $X$¶

Our probability distribution of $X$, shown as a table:

In [3]:
# our random variable X
dist_df = pd.DataFrame({"x": [3, 4, 6, 8],
                        "P(X = x)": [0.1, 0.2, 0.4, 0.3]})
dist_df
Out[3]:
x P(X = x)
0 3 0.1
1 4 0.2
2 6 0.4
3 8 0.3
In [4]:
plot_dist(dist_df, save=True)

Let's use this probability distribution to generate a table of $X(s)$, i.e., random variable values for many many samples.

In [5]:
# copied from above
# def simulate_samples(df, xname="x", pname="P(X = x)", size=1):
#     return np.random.choice(
#                 df[xname], # draw from these choiecs
#                 size=size, # this many times
#                 p=df[pname]) # according to this distribution

N = 80000
all_samples = simulate_samples(dist_df, size=N)
sim_df = pd.DataFrame({"X(s)": all_samples})
sim_df
Out[5]:
X(s)
0 8
1 6
2 6
3 8
4 3
... ...
79995 4
79996 6
79997 8
79998 4
79999 8

80000 rows × 1 columns



Let's check how well this simulated sample matches our probability distribution!

In [6]:
plot_simulated_dist(sim_df, "X(s)")
plt.title("Simulated distribution of $X$")
plt.show()
X(s)
E[•] 5.894812
Var(•) 2.867059
SD(•) 1.693239
In [7]:
# the tabular view of the above plot
sim_df.value_counts("X(s)").sort_values()/N
Out[7]:
X(s)
3    0.099138
4    0.200175
8    0.296287
6    0.404400
dtype: float64



Die Is the Singular; Dice Is the Plural¶

Let X be the outcome of a single die roll. X is a random variable.

In [4]:
# our random variable X
roll_df = pd.DataFrame({"x": [1, 2, 3, 4, 5, 6],
                        "P(X = x)": np.ones(6)/6})
roll_df
Out[4]:
x P(X = x)
0 1 0.166667
1 2 0.166667
2 3 0.166667
3 4 0.166667
4 5 0.166667
5 6 0.166667
In [5]:
plot_dist(roll_df)



Sum of 2 Dice Rolls¶

Here's the distribution of a single die roll:

In [9]:
roll_df = pd.DataFrame({"x": [1, 2, 3, 4, 5, 6],
                        "P(X = x)": np.ones(6)/6})
roll_df
Out[9]:
x P(X = x)
0 1 0.166667
1 2 0.166667
2 3 0.166667
3 4 0.166667
4 5 0.166667
5 6 0.166667

Let $X_1, X_2$ are the outcomes of two dice rolls. Note $X_1$ and $X_2$ are i.i.d. (independent and identically distributed).

Below I call a helper function simulate_iid_df, which simulates an 80,000-row table of $X_1, X_2$ values. It uses np.random.choice(arr, size, p) link where arr is the array the values and p is the probability associated with choosing each value. If you're interested in the implementation details, scroll up.

In [8]:
N = 80000
sim_rolls_df = simulate_iid_df(roll_df, nvars=2, rows=N)
sim_rolls_df
Out[8]:
X_1 X_2
0 2 4
1 5 4
2 3 4
3 6 6
4 2 5
... ... ...
79995 6 5
79996 1 4
79997 4 3
79998 6 4
79999 5 1

80000 rows × 2 columns

Define the following random variables, which are functions of $X_1$ and $X_2$:

  • $Y = X_1 + X_1 = 2 X_1$
  • $Z = X_1 + X_2$

We can use our simulated values of $X_1, X_2$ to create new columns $Y$ and $Z$:

In [12]:
sim_rolls_df['Y'] = 2 * sim_rolls_df['X_1']
sim_rolls_df['Z'] = sim_rolls_df['X_1'] + sim_rolls_df['X_2']
sim_rolls_df
Out[12]:
X_1 X_2 Y Z
0 2 4 4 6
1 5 4 10 9
2 3 4 6 7
3 6 6 12 12
4 2 5 4 7
... ... ... ... ...
79995 6 5 12 11
79996 1 4 2 5
79997 4 3 8 7
79998 6 4 12 10
79999 5 1 10 6

80000 rows × 4 columns

Now that we have simulated samples of $Y$ and $Z$, we can plot histograms to see their distributions!


Distribution of $Y$, which was twice the value of our first die roll:

In [13]:
plot_simulated_dist(sim_rolls_df, "Y", save=True)
Y
E[•] 6.999625
Var(•) 11.688446
SD(•) 3.418837

Distribution of $Z$, the sum of two IID dice rolls:

In [14]:
plot_simulated_dist(sim_rolls_df, "Z", save=True, color='gold')
Z
E[•] 7.003800
Var(•) 5.848734
SD(•) 2.418416

Let's compare the expectations and variances of these simulated distributions of $Y$ and $Z$.

  • We computed:
    • $\mathbb{E}[Y]$ as np.mean(sim_rolls_df['Y'])
    • $\text{Var}(Y)$ as np.var(sim_rolls_df['Y']
    • etc.
  • The larger your simulated rows $N$, the closer the simulated expectation will be to the true expectation.
  • Our approach is tedious--we have to simulate an entire population, then reduce it down to expectation/variance/standard deviation. There has to be a better way!
In [27]:
stats_df_multi(sim_rolls_df, ["Y", "Z"])
Out[27]:
Y Z
E[•] 7.013375 7.002725
Var(•) 11.722368 5.870866
SD(•) 3.423794 2.422987



Which would you pick?¶

  • $\large Y_A = 10 X_1 + 10 X_2 $
  • $\large Y_B = \sum\limits_{i=1}^{20} X_i$
  • $\large Y_C = 20 X_1$

First let's construct the probability distribution for a single coin. This will let us flip 20 IID coins later.

In [16]:
# First construct probability distribution for a single fair coin
p = 0.5
coin_df = pd.DataFrame({"x": [0, 1], # [Tails, Heads]
                        "P(X = x)": [p, 1 - p]})
coin_df
Out[16]:
x P(X = x)
0 0 0.5
1 1 0.5

Choice A:¶

$\large Y_A = 10 X_1 + 10 X_2 $

In [36]:
# Flip 20 iid coins, each exactly once
flips20_df = simulate_iid_df(coin_df, nvars=20, rows=1)

# Construct Y_A from this sample
flips20_df["Y_A"] = 10*flips20_df["X_1"] + 10*flips20_df["X_2"]

print("Sample of size 1:")
display(flips20_df)
print("Y_A:", flips20_df.loc[0,"Y_A"])

Sample of size 1:
X_1 X_2 X_3 X_4 X_5 X_6 X_7 X_8 X_9 X_10 ... X_12 X_13 X_14 X_15 X_16 X_17 X_18 X_19 X_20 Y_A
0 0 0 0 0 1 0 0 1 1 1 ... 0 0 1 1 1 0 1 0 1 0

1 rows × 21 columns

Y_A: 0

Choice B:¶

$\large Y_B = \sum\limits_{i=1}^{20} X_i$

In [37]:
# Flip 20 iid coins, each exactly once
flips20_df = simulate_iid_df(coin_df, nvars=20, rows=1)

# Construct Y_B from this sample
flips20_df["Y_B"] = flips20_df.sum(axis=1) # sum all coins

display(flips20_df)
print("Y_B:", flips20_df.loc[0,"Y_B"])
X_1 X_2 X_3 X_4 X_5 X_6 X_7 X_8 X_9 X_10 ... X_12 X_13 X_14 X_15 X_16 X_17 X_18 X_19 X_20 Y_B
0 0 1 0 1 0 0 1 1 0 0 ... 1 0 0 1 1 0 0 0 0 8

1 rows × 21 columns

Y_B: 8

Choice C:¶

$\large Y_C = 20 X_1$

In [43]:
# Flip 20 iid coins, each exactly once
flips20_df = simulate_iid_df(coin_df, nvars=20, rows=1)

# Construct Y_C from this sample
flips20_df["Y_C"] = 20*flips20_df["X_1"]

display(flips20_df[["X_1", "Y_C"]])
print("Y_C:", flips20_df.loc[0,"Y_C"])
X_1 Y_C
0 0 0 
Y_C: 0


If you're curious as to what these distributions look like, I've simulated some populations:

In [35]:
# simulate one big population of 20 fair flips, N = 80,000
N = 80000
df_coins = simulate_iid_df(coin_df, nvars=20, rows=N)
print(f"{N} simulated samples")


# construct Y_A, Y_B, Y_C from this population
df_coins["Y_A"] = 10*df_coins["X_1"] + 10*df_coins["X_2"]
df_coins["Y_B"] = df_coins.loc[:,"X_1":"X_20"].sum(axis=1)
df_coins["Y_C"] = 5 + 15*df_coins["X_1"]
plot_simulated_dist_multi(df_coins, ["Y_A", "Y_B", "Y_C"])

# adjust axes for nicer plotting
axs = plt.gcf().axes
axs[1].set_xlim(axs[0].get_xlim())    # Y_B
axs[1].set_xticks([0, 5, 10, 15, 20]) # Y_B
axs[2].set_xlim(axs[0].get_xlim())    # Y_C
plt.show()

80000 simulated samples
Y_A Y_B Y_C
E[•] 10.012250 10.014463 12.524563
Var(•) 50.015475 4.990516 56.250100
SD(•) 7.072162 2.233946 7.500007



If we flipped 100 coins, look how beautiful the Binomial distribution looks:

In [33]:
# simulate one big population of 100 fair flips, N = 80,000
N = 80000
df_coins = simulate_iid_df(coin_df, nvars=50, rows=N)
df_coins["Y_100"] = df_coins.loc[:,"X_1":"X_50"].sum(axis=1)
plot_simulated_dist(df_coins, "Y_100")
plt.xticks([0, 10, 20, 30, 40, 50])
plt.show()
Y_100
E[•] 25.005425
Var(•) 12.545202
SD(•) 3.541921




From Population to Sample¶

Remember the population distribution we looked at earlier:

In [61]:
dist_df
Out[61]:
x P(X = x)
0 3 0.1
1 4 0.2
2 6 0.4
3 8 0.3
In [96]:
# a population generated from the distribution
N = 100000
all_samples = simulate_samples(dist_df, size=N)
sim_pop_df = pd.DataFrame({"X(s)": all_samples})
sim_pop_df
Out[96]:
X(s)
0 6
1 8
2 8
3 4
4 8
... ...
99995 6
99996 8
99997 8
99998 8
99999 8

100000 rows × 1 columns




Suppose we draw a sample of size 100 from this giant population.

We are performing Random Sampling with Replacement: df.sample(n, replace=True) (link)

In [84]:
n = 100      # size of our sample
sample_df = (
             sim_pop_df.sample(n, replace=True)
             
             # some reformatting below
             .reset_index(drop=True)
             .rename(columns={"X(s)": "X"})
            )
sample_df
Out[84]:
X
0 6
1 8
2 6
3 6
4 3
... ...
95 8
96 6
97 6
98 3
99 8

100 rows × 1 columns

Our sample distribution (n = 100):

In [101]:
sns.histplot(data=sample_df, x='X', stat='probability')
plt.xticks([3, 4, 6, 8])
plt.title(f"Sample (n = 100)")
plt.show()

print("Mean of Sample:", np.mean(sample_df['X']))

Mean of Sample: 5.77


Compare this to our original population (N = 80,000):

In [97]:
plot_simulated_dist(sim_df, "X(s)")
plt.title("Population of $X$")
plt.show()
X(s)
E[•] 5.893600
Var(•) 2.888465
SD(•) 1.699549