ch245 — Conditional Probability

import numpy as np

# Uniform distribution on d6
Omega = set(range(1, 7))

def P(event, omega=Omega):
    return len(event) / len(omega)

def P_given(A, B, omega=Omega):
    """P(A | B) = P(A ∩ B) / P(B)"""
    if len(B) == 0:
        raise ValueError("Cannot condition on impossible event P(B)=0")
    return len(A & B) / len(B)  # equivalent for uniform: |A∩B|/|B|

E_even = {2, 4, 6}
E_six = {6}
E_prime = {2, 3, 5}

print(f"P(6) = {P(E_six):.4f}")
print(f"P(6 | even) = {P_given(E_six, E_even):.4f}")
print(f"P(prime | even) = {P_given(E_prime, E_even):.4f}")
print(f"P(even | prime) = {P_given(E_even, E_prime):.4f}")

# Geometric view: conditioning restricts the sample space
import matplotlib.pyplot as plt
import matplotlib.patches as patches

fig, axes = plt.subplots(1, 2, figsize=(10, 4))

for ax, highlight_B, title in zip(
    axes,
    [False, True],
    ["Full Ω — P(A)", "Restricted to B — P(A|B)"]
):
    ax.set_xlim(-0.5, 6.5)
    ax.set_ylim(-0.5, 1.5)
    ax.set_aspect('equal')
    ax.set_title(title, fontsize=11)
    ax.axis('off')

    # Draw all outcomes as circles
    for i, val in enumerate(range(1, 7), start=0):
        is_even = val % 2 == 0
        is_six = val == 6
        
        if highlight_B and not is_even:
            color = 'lightgray'  # outside B — grayed out
        elif is_six:
            color = 'tomato'     # A ∩ B
        elif is_even:
            color = 'steelblue'  # B but not A
        else:
            color = 'white'
        
        circle = plt.Circle((i, 0.5), 0.35, color=color, ec='black', lw=1.5)
        ax.add_patch(circle)
        ax.text(i, 0.5, str(val), ha='center', va='center', fontsize=14, fontweight='bold')

    ax.text(3, -0.2, 'Red = outcome 6   Blue = other evens   Gray = excluded', 
            ha='center', fontsize=8, color='gray')

plt.tight_layout()
plt.show()

# Draw two aces from a standard 52-card deck without replacement
# P(1st is ace) = 4/52
# P(2nd is ace | 1st was ace) = 3/51
# P(both aces) = P(1st ace) * P(2nd ace | 1st ace)

p_first_ace = 4 / 52
p_second_ace_given_first = 3 / 51
p_both_aces = p_first_ace * p_second_ace_given_first

print(f"P(first card is ace) = {p_first_ace:.6f}")
print(f"P(second is ace | first was ace) = {p_second_ace_given_first:.6f}")
print(f"P(both aces) = {p_both_aces:.6f}")
print(f"  = 1 in {1/p_both_aces:.1f} deals")

# Simulation
rng = np.random.default_rng(seed=42)
n_trials = 1_000_000
count_both_aces = 0

deck = np.array([0]*4 + [1]*48)  # 0=ace, 1=non-ace
for _ in range(n_trials):
    draw = rng.choice(deck, size=2, replace=False)
    if draw[0] == 0 and draw[1] == 0:
        count_both_aces += 1

print(f"\nSimulation ({n_trials:,} trials): {count_both_aces/n_trials:.6f}")

# Factory example: two machines produce parts
# Machine 1 produces 60% of parts, has 2% defect rate
# Machine 2 produces 40% of parts, has 5% defect rate
# What is the overall defect rate?

p_B1 = 0.60  # P(from machine 1)
p_B2 = 0.40  # P(from machine 2)
p_defect_given_B1 = 0.02
p_defect_given_B2 = 0.05

# Law of total probability
p_defect = p_defect_given_B1 * p_B1 + p_defect_given_B2 * p_B2

print("Law of Total Probability: P(defect) = P(defect|B1)P(B1) + P(defect|B2)P(B2)")
print(f"  = {p_defect_given_B1} × {p_B1} + {p_defect_given_B2} × {p_B2}")
print(f"  = {p_defect_given_B1 * p_B1:.4f} + {p_defect_given_B2 * p_B2:.4f}")
print(f"  = {p_defect:.4f}  ({p_defect*100:.2f}% overall defect rate)")

# Verify: two dice rolls are independent
from itertools import product

Omega2 = list(product(range(1, 7), range(1, 7)))

def P2_given(A, B):
    """P(A|B) for uniform distribution over Omega2."""
    A, B = set(A), set(B)
    return len(A & B) / len(B)

E_d1_six = {o for o in Omega2 if o[0] == 6}
E_d2_six = {o for o in Omega2 if o[1] == 6}

p_d1_six = len(E_d1_six) / len(Omega2)
p_d1_six_given_d2_six = P2_given(E_d1_six, E_d2_six)

print(f"P(d1=6) = {p_d1_six:.4f}")
print(f"P(d1=6 | d2=6) = {p_d1_six_given_d2_six:.4f}")
print(f"Equal → independent: {abs(p_d1_six - p_d1_six_given_d2_six) < 1e-10}")