# --- Setup: load our helpers and the data --------------------------------
import sys, pathlib
for _p in [pathlib.Path.cwd(), pathlib.Path.cwd().parent]:
    if (_p / 'lib' / 'statslab.py').exists():
        sys.path.insert(0, str(_p / 'lib')); break
import statslab as sl

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
sl.use_course_style()
rng = np.random.default_rng(12)  # fixed seed -> reproducible; learner can change it
print("Ready.")

Ready.

faces = np.array([1, 2, 3, 4, 5, 6])
probs = np.full(6, 1/6)          # each face has probability 1/6

print("Outcome  Probability")
for x, p in zip(faces, probs):
    print(f"   {x}        {p:.4f}")
print(f"\nProbabilities sum to {probs.sum():.4f}  (must be 1).")

fig, ax = plt.subplots()
ax.bar(faces, probs, color='#4c72b0', edgecolor='white')
ax.set_xlabel('die face (value of X)'); ax.set_ylabel('probability P(X = x)')
ax.set_title('Distribution of a fair die — six equal bars')
ax.set_ylim(0, 0.25)
plt.show()

Outcome  Probability
   1        0.1667
   2        0.1667
   3        0.1667
   4        0.1667
   5        0.1667
   6        0.1667

Probabilities sum to 1.0000  (must be 1).

E = np.sum(faces * probs)        # sum of x * P(x)
print(f"E[X] = {E}   (the textbook value is 3.5)")

E[X] = 3.5   (the textbook value is 3.5)

Var = np.sum((faces - E)**2 * probs)   # sum of (x - mu)^2 * P(x)
SD  = np.sqrt(Var)
print(f"Var[X] = {Var:.4f}   (textbook 2.9167)")
print(f"SD[X]  = {SD:.4f}   (textbook 1.7078)")

Var[X] = 2.9167   (textbook 2.9167)
SD[X]  = 1.7078   (textbook 1.7078)

n_rolls = 5000                              # <- tinker: try 200, then 50000
rolls = rng.integers(1, 7, size=n_rolls)    # whole numbers 1..6, each equally likely

running_avg = np.cumsum(rolls) / np.arange(1, n_rolls + 1)

fig, ax = plt.subplots()
ax.plot(np.arange(1, n_rolls + 1), running_avg, color='#4c72b0', lw=1.2,
        label='running average of rolls')
ax.axhline(E, color='black', lw=2, ls='--', label=f'E[X] = {E}')
ax.set_xscale('log')                        # log x-axis spreads out the early, jumpy part
ax.set_xlabel('number of rolls so far (log scale)'); ax.set_ylabel('average value rolled')
ax.set_title('Law of large numbers: the sample mean converges to E[X]')
ax.legend(); plt.show()

print(f"Average of first 10 rolls:    {running_avg[9]:.3f}")
print(f"Average of all {n_rolls} rolls: {running_avg[-1]:.3f}   (target: {E})")

Average of first 10 rolls:    3.200
Average of all 5000 rolls: 3.509   (target: 3.5)

sim_var = rolls.var(ddof=0)        # population-style variance of the simulated rolls
print(f"Simulated variance of {n_rolls} rolls: {sim_var:.4f}")
print(f"Exact Var[X] for a fair die:        {Var:.4f}")

Simulated variance of 5000 rolls: 2.9387
Exact Var[X] for a fair die:        2.9167

win_amt, lose_amt = 10, -2
p_win = 1/6
E_payoff = p_win * win_amt + (1 - p_win) * lose_amt
E_payoff = round(E_payoff, 10) + 0.0   # round off a tiny float residual so 0 prints as 0.00
print(f"E[payoff] = (1/6)*{win_amt} + (5/6)*({lose_amt}) = ${E_payoff:.2f}  -> a fair game")

# Simulate playing the game many times and average the result.
n_plays = 100_000                       # <- tinker: try a small number like 20
die = rng.integers(1, 7, size=n_plays)
payoffs = np.where(die == 6, win_amt, lose_amt)
print(f"Average payoff over {n_plays:,} plays: ${payoffs.mean():.3f}  (heading to $0.00)")

E[payoff] = (1/6)*10 + (5/6)*(-2) = $0.00  -> a fair game
Average payoff over 100,000 plays: $-0.019  (heading to $0.00)

births = sl.load_csv('births14.csv', url='https://www.openintro.org/data/csv/births14.csv')
print(f"Rows before dropping missing visits: {len(births)}")

visits = births['visits'].dropna()
n = len(visits)
print(f"Rows after dropna:                   {n}")
print(f"mean   = {visits.mean():.2f}")
print(f"var (sample, ddof=1) = {visits.var(ddof=1):.2f}")
print(f"min = {visits.min():.0f}   max = {visits.max():.0f}")

Rows before dropping missing visits: 1000
Rows after dropna:                   944
mean   = 11.35
var (sample, ddof=1) = 16.88
min = 0   max = 30

counts = visits.value_counts().sort_index()      # how many births at each k
k_vals = counts.index.to_numpy(dtype=int)        # the values 0,1,2,...,30
P = counts.to_numpy() / n                         # P(X = k)

dist = pd.DataFrame({'k (visits)': k_vals, 'count': counts.to_numpy(), 'P(X=k)': P})
print(f"Probabilities sum to {P.sum():.4f}  (must be 1).")
print("\nA few rows of the empirical distribution:")
print(dist.head(12).to_string(index=False))
print("...")
print(dist.tail(4).to_string(index=False))

Probabilities sum to 1.0000  (must be 1).

A few rows of the empirical distribution:
 k (visits)  count   P(X=k)
          0     17 0.018008
          1      6 0.006356
          2      4 0.004237
          3     13 0.013771
          4     14 0.014831
          5     15 0.015890
          6     27 0.028602
          7     37 0.039195
          8     51 0.054025
          9     56 0.059322
         10    132 0.139831
         11     86 0.091102
...
 k (visits)  count   P(X=k)
         27      1 0.001059
         28      1 0.001059
         29      1 0.001059
         30      2 0.002119

fig, ax = plt.subplots()
ax.bar(k_vals, P, color='#4c72b0', edgecolor='white')
ax.set_xlabel('number of prenatal visits (k)'); ax.set_ylabel('probability P(X = k)')
ax.set_title('Empirical distribution of prenatal visits (n = %d)' % n)
plt.show()

E_visits   = np.sum(k_vals * P)                       # sum k * P(k)
Var_visits = np.sum((k_vals - E_visits)**2 * P)       # sum (k - E)^2 * P(k)
SD_visits  = np.sqrt(Var_visits)

print(f"E[X] from distribution = {E_visits:.4f}")
print(f"data mean of visits    = {visits.mean():.4f}   <- should match E[X]")
print()
print(f"Var[X] from distribution      = {Var_visits:.4f}")
print(f"visits.var(ddof=0) (population) = {visits.var(ddof=0):.4f}   <- should match Var[X]")
print(f"visits.var(ddof=1) (sample)     = {visits.var(ddof=1):.4f}   <- a hair larger")
print()
print(f"SD[X] = {SD_visits:.4f} visits")

E[X] from distribution = 11.3517
data mean of visits    = 11.3517   <- should match E[X]

Var[X] from distribution      = 16.8594
visits.var(ddof=0) (population) = 16.8594   <- should match Var[X]
visits.var(ddof=1) (sample)     = 16.8772   <- a hair larger

SD[X] = 4.1060 visits

n_draws = 4000                                   # <- tinker: try 100, then 20000
sample = rng.choice(visits.to_numpy(), size=n_draws, replace=True)
running = np.cumsum(sample) / np.arange(1, n_draws + 1)

fig, ax = plt.subplots()
ax.plot(np.arange(1, n_draws + 1), running, color='#009E73', lw=1.2,
        label='running average of sampled births')
ax.axhline(E_visits, color='black', lw=2, ls='--', label=f'E[X] = {E_visits:.2f}')
ax.set_xscale('log')
ax.set_xlabel('number of births drawn so far (log scale)'); ax.set_ylabel('average visits')
ax.set_title('Law of large numbers on real data: sample mean -> E[X]')
ax.legend(); plt.show()

print(f"Average of first 10 draws:   {running[9]:.2f}")
print(f"Average of all {n_draws} draws: {running[-1]:.2f}   (target E[X] = {E_visits:.2f})")

Average of first 10 draws:   11.30
Average of all 4000 draws: 11.31   (target E[X] = 11.35)

Lesson 7 — Random variables, expectation & variance¶

1) What is a random variable?¶

2) A fair six-sided die — a distribution we know exactly¶

3) Expected value E[X] — the "center of gravity"¶

4) Variance Var[X] and standard deviation — how spread out is it?¶

5) The law of large numbers — watch the average converge¶

6) Expected value of a bet — is the game fair?¶

7) A REAL random variable — number of prenatal visits¶

8) Build the empirical distribution: P(X = k) = count of k ÷ n¶

9) Compute E[X] and Var[X] FROM the distribution¶

10) Law of large numbers again — but on real data¶

Now you try¶

What you learned¶