import numpy as np, pandas as pd, matplotlib.pyplot as plt
from itertools import product
rng = np.random.default_rng(30)
BLUE="#2563eb"; INDIGO="#4f46e5"; CYAN="#0891b2"; AMBER="#d97706"; GREEN="#059669"; PURPLE="#7c3aed"; PINK="#db2777"; INK="#1a2138"; GRID="#e6e9f2"
plt.rcParams.update({"figure.facecolor":"white","axes.facecolor":"white","figure.dpi":110,"font.size":11,
   "axes.edgecolor":GRID,"axes.grid":True,"grid.color":GRID,"axes.axisbelow":True,"axes.spines.top":False,
   "axes.spines.right":False,"axes.titlesize":12,"axes.titleweight":"bold","legend.frameon":False})
print("ready")

ready

sample_space = list(product(range(1,7), range(1,7)))
print("rolling two dice -> sample space size:", len(sample_space))
print("first few outcomes:", sample_space[:5])

# events are subsets of the sample space
sum7      = [o for o in sample_space if sum(o)==7]
doubles   = [o for o in sample_space if o[0]==o[1]]
sum_gt9   = [o for o in sample_space if sum(o)>9]
print("\nevent {sum = 7}   :", len(sum7), "outcomes ->", sum7)
print("event {doubles}   :", len(doubles), "outcomes")
print("event {sum > 9}   :", len(sum_gt9), "outcomes")

rolling two dice -> sample space size: 36
first few outcomes: [(1, 1), (1, 2), (1, 3), (1, 4), (1, 5)]

event {sum = 7}   : 6 outcomes -> [(1, 6), (2, 5), (3, 4), (4, 3), (5, 2), (6, 1)]
event {doubles}   : 6 outcomes
event {sum > 9}   : 6 outcomes

totals = [sum(o) for o in sample_space]
dist = pd.Series(totals).value_counts().sort_index()
prob = (dist/36).round(3)
print("P(sum = k) for two dice:")
for k in dist.index: print(f"   {k:>2}: {dist[k]}/36 = {prob[k]:.3f}")

fig, ax = plt.subplots(figsize=(7.5,3.6))
ax.bar(dist.index, dist.values/36, color=BLUE, edgecolor="white")
ax.set_xticks(range(2,13)); ax.set_xlabel("sum of two dice"); ax.set_ylabel("probability")
ax.set_title("classical probability: P(sum) = favorable / 36"); plt.tight_layout(); plt.show()
print("\nmost likely sum is 7 (P = 6/36 = 0.167); 2 and 12 are rarest (1/36 each).")

P(sum = k) for two dice:
    2: 1/36 = 0.028
    3: 2/36 = 0.056
    4: 3/36 = 0.083
    5: 4/36 = 0.111
    6: 5/36 = 0.139
    7: 6/36 = 0.167
    8: 5/36 = 0.139
    9: 4/36 = 0.111
   10: 3/36 = 0.083
   11: 2/36 = 0.056
   12: 1/36 = 0.028

most likely sum is 7 (P = 6/36 = 0.167); 2 and 12 are rarest (1/36 each).

N = 4000
flips = rng.integers(0, 2, N)            # 1 = heads
run_heads = np.cumsum(flips) / np.arange(1, N+1)
rolls = rng.integers(1, 7, N)
run_six   = np.cumsum(rolls==6) / np.arange(1, N+1)

fig, ax = plt.subplots(1, 2, figsize=(11,3.8))
ax[0].plot(run_heads, color=GREEN); ax[0].axhline(0.5, ls="--", color=PINK)
ax[0].set_title(f"coin: proportion of heads (-> 0.5)"); ax[0].set_xlabel("number of flips"); ax[0].set_ylim(0,1)
ax[1].plot(run_six, color=INDIGO); ax[1].axhline(1/6, ls="--", color=PINK)
ax[1].set_title("die: proportion of sixes (-> 1/6)"); ax[1].set_xlabel("number of rolls"); ax[1].set_ylim(0,0.5)
plt.tight_layout(); plt.show()
print(f"after {N} flips: P(heads) approx {run_heads[-1]:.3f}  (theory 0.500)")
print(f"after {N} rolls: P(six)   approx {run_six[-1]:.3f}  (theory {1/6:.3f})")

after 4000 flips: P(heads) approx 0.509  (theory 0.500)
after 4000 rolls: P(six)   approx 0.167  (theory 0.167)

# P(at least one six in two rolls): the hard way vs the complement
p_no_six_one_roll = 5/6
p_at_least_one = 1 - p_no_six_one_roll**2
print(f"P(no six in two rolls)      = (5/6)^2 = {p_no_six_one_roll**2:.3f}")
print(f"P(at least one six)         = 1 - that = {p_at_least_one:.3f}")

# verify empirically
trials = rng.integers(1,7,(100_000,2))
emp = np.mean((trials==6).any(axis=1))
print(f"empirical (100k trials)     = {emp:.3f}")

events = {"impossible (sum=13)":0.0, "two sixes":1/36, "a seven":6/36,
          "not a seven":30/36, "certain (sum<=12)":1.0}
fig, ax = plt.subplots(figsize=(8,1.7))
for name,p in events.items(): ax.scatter(p, 0, s=80, zorder=3)
for name,p in events.items(): ax.annotate(name, (p,0), rotation=30, fontsize=7.5, ha="left", va="bottom", xytext=(0,6), textcoords="offset points")
ax.axhline(0, color=INK, lw=2); ax.set_xlim(-0.05,1.15); ax.set_yticks([]); ax.set_xticks([0,0.5,1])
ax.set_xlabel("probability"); ax.set_title("the probability scale: 0 = impossible, 1 = certain"); plt.tight_layout(); plt.show()

P(no six in two rolls)      = (5/6)^2 = 0.694
P(at least one six)         = 1 - that = 0.306

empirical (100k trials)     = 0.307

flips = rng.integers(0,2,200_000)        # 1 = heads
# find positions right after a run of >=3 consecutive heads
after_streak = []
streak = 0
for i in range(len(flips)-1):
    streak = streak+1 if flips[i]==1 else 0
    if streak >= 3:
        after_streak.append(flips[i+1])
after_streak = np.array(after_streak)
print(f"flips following a run of >=3 heads: {len(after_streak):,}")
print(f"proportion that were ALSO heads  : {after_streak.mean():.3f}")
print("theory: still 0.500 -- the coin does not remember the streak.")

flips following a run of >=3 heads: 24,950
proportion that were ALSO heads  : 0.500
theory: still 0.500 -- the coin does not remember the streak.

⚙️ Setup¶