Random Variables and Distributions

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102

import numpy as np
from scipy import stats
 
def demonstrate_cdf_properties():
    """
    Demonstrate the fundamental properties of CDFs.
    """
    print("CDF Properties Demonstration")
    print("=" * 60)
    
    # Example 1: Discrete CDF (Binomial)
    print("
1. Discrete CDF: Binomial(n=10, p=0.5)")
    print("-" * 60)
    
    binom = stats.binom(n=10, p=0.5)
    
    # Show step function nature
    print("
Step function behavior:")
    for x in [2.0, 2.5, 2.9, 3.0, 3.1, 3.5, 4.0]:
        print(f"  F({x:.1f}) = P(X ≤ {x:.1f}) = {binom.cdf(x):.6f}")
    
    print("
Note: F is constant between integers (steps only at support)")
    
    # Example 2: Continuous CDF (Normal)
    print("
2. Continuous CDF: Normal(μ=0, σ=1)")
    print("-" * 60)
    
    normal = stats.norm(0, 1)
    
    # Smooth behavior
    print("
Smooth behavior (no jumps):")
    for x in np.linspace(-2, 2, 9):
        print(f"  F({x:+5.2f}) = {normal.cdf(x):.6f}")
    
    # Example 3: Verify monotonicity
    print("
3. Monotonicity Check")
    print("-" * 60)
    
    x_values = np.linspace(-3, 3, 100)
    cdf_values = normal.cdf(x_values)
    
    # Check that each value >= previous
    is_monotonic = np.all(np.diff(cdf_values) >= 0)
    print(f"CDF is monotonically non-decreasing: {is_monotonic}")
    
    # Example 4: Boundary conditions
    print("
4. Boundary Conditions")
    print("-" * 60)
    
    print(f"lim(x → -∞) F(x) = {normal.cdf(-100):.10f} ≈ 0")
    print(f"lim(x → +∞) F(x) = {normal.cdf(100):.10f} ≈ 1")
 
 
def cdf_for_probability_computations():
    """
    Show how CDFs simplify probability computations.
    """
    print("
 
Using CDFs for Probability Computation")
    print("=" * 60)
    
    normal = stats.norm(0, 1)
    
    print("
Key formulas:")
    print("  P(X ≤ a) = F(a)")
    print("  P(X > a) = 1 - F(a)")
    print("  P(a < X ≤ b) = F(b) - F(a)")
    print("  P(a ≤ X ≤ b) = F(b) - F(a) [for continuous X]")
    
    a, b = -1, 2
    print(f"
Example with N(0,1), a={a}, b={b}:")
    print(f"  P(X ≤ {a}) = F({a}) = {normal.cdf(a):.6f}")
    print(f"  P(X > {b}) = 1 - F({b}) = {1 - normal.cdf(b):.6f}")
    print(f"  P({a} < X ≤ {b}) = F({b}) - F({a}) = {normal.cdf(b) - normal.cdf(a):.6f}")
    
    # For discrete: P(X = k) = F(k) - F(k-1)
    print("
For discrete RVs:")
    print("  P(X = k) = F(k) - F(k-1)") 
    
    binom = stats.binom(10, 0.5)
    k = 5
    p_k_via_pmf = binom.pmf(k)
    p_k_via_cdf = binom.cdf(k) - binom.cdf(k-1)
    print(f"
Binomial(10, 0.5): P(X = {k})")
    print(f"  Via PMF: {p_k_via_pmf:.6f}")
    print(f"  Via CDF diff: F({k}) - F({k-1}) = {p_k_via_cdf:.6f}")
 
 
demonstrate_cdf_properties()
cdf_for_probability_computations()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107

import numpy as np
from scipy import stats, integrate
 
def discrete_relationships():
    """
    Demonstrate PMF ↔ CDF relationships for discrete RVs.
    """
    print("Discrete RV: PMF ↔ CDF Relationships")
    print("=" * 60)
    
    # Poisson distribution
    lam = 3
    poisson = stats.poisson(lam)
    
    print(f"
Poisson(λ={lam})")
    print("
Computing CDF from PMF by summation:")
    print(f"{'k':>4} {'p(k)':>10} {'F(k) from sum':>15} {'F(k) from scipy':>18}")
    print("-" * 50)
    
    cumulative = 0
    for k in range(10):
        pmf_k = poisson.pmf(k)
        cumulative += pmf_k
        cdf_k = poisson.cdf(k)
        print(f"{k:>4} {pmf_k:>10.6f} {cumulative:>15.6f} {cdf_k:>18.6f}")
    
    print("
Recovering PMF from CDF:")
    print("  p(k) = F(k) - F(k-1)")
    
    for k in range(1, 6):
        recovered_pmf = poisson.cdf(k) - poisson.cdf(k-1)
        actual_pmf = poisson.pmf(k)
        print(f"  p({k}) = F({k}) - F({k-1}) = {recovered_pmf:.6f} (actual: {actual_pmf:.6f})")
 
 
def continuous_relationships():
    """
    Demonstrate PDF ↔ CDF relationships for continuous RVs.
    """
    print("
 
Continuous RV: PDF ↔ CDF Relationships")
    print("=" * 60)
    
    # Standard Normal
    normal = stats.norm(0, 1)
    
    print("
Normal(0, 1)")
    print("
Computing CDF from PDF by integration:")
    print(f"{'x':>6} {'f(x)':>10} {'∫f(t)dt from -∞':>20} {'F(x) from scipy':>18}")
    print("-" * 56)
    
    for x in [-2, -1, 0, 1, 2]:
        pdf_x = normal.pdf(x)
        # Numerical integration
        cdf_integral, _ = integrate.quad(normal.pdf, -np.inf, x)
        cdf_scipy = normal.cdf(x)
        print(f"{x:>6.1f} {pdf_x:>10.6f} {cdf_integral:>20.10f} {cdf_scipy:>18.10f}")
    
    print("
Recovering PDF from CDF by numerical differentiation:")
    print("  f(x) ≈ [F(x+h) - F(x-h)] / (2h)")
    
    h = 0.0001
    for x in [-1.0, 0.0, 1.0, 2.0]:
        numerical_derivative = (normal.cdf(x + h) - normal.cdf(x - h)) / (2 * h)
        actual_pdf = normal.pdf(x)
        print(f"  f({x:+.1f}) ≈ {numerical_derivative:.6f} (actual: {actual_pdf:.6f})")
 
 
def survival_function():
    """
    Introduce the survival function (complement of CDF).
    """
    print("
 
The Survival Function: S(x) = 1 - F(x) = P(X > x)")
    print("=" * 60)
    
    print("
The survival function is critical in:")
    print("  - Reliability engineering: P(component survives past time t)")
    print("  - Medical studies: P(patient survives past time t)")
    print("  - ML: P-values, tail probabilities, outlier detection")
    
    normal = stats.norm(0, 1)
    
    print("
For N(0, 1):")
    for x in [1, 1.645, 1.96, 2.576, 3]:
        survival = normal.sf(x)  # sf = survival function
        complement = 1 - normal.cdf(x)
        print(f"  S({x:.3f}) = P(X > {x:.3f}) = {survival:.6f}")
    
    print("
Note: scipy provides both cdf(x) and sf(x) = 1 - cdf(x)")
    print("      sf(x) is more numerically stable for small tail probs")
 
 
discrete_relationships()
continuous_relationships()
survival_function()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140

import numpy as np
from scipy import stats
 
def demonstrate_quantiles():
    """
    Demonstrate quantile functions and inverse transform sampling.
    """
    print("Quantile Functions (Inverse CDF)")
    print("=" * 60)
    
    # Standard Normal quantiles
    normal = stats.norm(0, 1)
    
    print("
Standard Normal N(0, 1) Quantiles:")
    print("-" * 40)
    
    quantiles = [0.001, 0.01, 0.025, 0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 0.95, 0.975, 0.99]
    for p in quantiles:
        q = normal.ppf(p)  # ppf = percent point function = inverse CDF
        print(f"  Q({p:.3f}) = {q:+.4f}")
    
    # Verify: F(Q(p)) = p
    print("
Verification: F(Q(p)) = p")
    for p in [0.25, 0.5, 0.75]:
        q = normal.ppf(p)
        recovered_p = normal.cdf(q)
        print(f"  F(Q({p})) = F({q:.4f}) = {recovered_p:.6f}")
 
 
def inverse_transform_sampling():
    """
    Demonstrate inverse transform sampling.
    """
    print("
 
Inverse Transform Sampling")
    print("=" * 60)
    print("
Method: To sample from distribution F,")
    print("  1. Sample U ~ Uniform(0, 1)")
    print("  2. Return Q(U) = F⁻¹(U)")
    
    np.random.seed(42)
    n_samples = 100000
    
    # Target: Exponential(λ=2)
    true_lambda = 2.0
    exponential = stats.expon(scale=1/true_lambda)
    
    # Method 1: Direct sampling (scipy)
    direct_samples = exponential.rvs(n_samples)
    
    # Method 2: Inverse transform
    uniform_samples = np.random.uniform(0, 1, n_samples)
    inverse_transform_samples = exponential.ppf(uniform_samples)
    
    print(f"
Target: Exponential(λ={true_lambda})")
    print(f"True mean = 1/λ = {1/true_lambda}")
    print(f"
Direct sampling mean: {direct_samples.mean():.4f}")
    print(f"Inverse transform mean: {inverse_transform_samples.mean():.4f}")
    
    print(f"
True variance = 1/λ² = {1/true_lambda**2}")
    print(f"Direct sampling var: {direct_samples.var():.4f}")
    print(f"Inverse transform var: {inverse_transform_samples.var():.4f}")
    
    # For Exponential, we can derive the inverse CDF analytically
    print("
 
Analytical Inverse Transform for Exponential:")
    print("  CDF: F(x) = 1 - exp(-λx)")
    print("  Solving F(x) = u for x:")
    print("  u = 1 - exp(-λx)")
    print("  exp(-λx) = 1 - u")
    print("  x = -log(1-u)/λ = Q(u)")
    
    # Verify
    analytical_samples = -np.log(1 - uniform_samples) / true_lambda
    print(f"
Analytical inverse transform mean: {analytical_samples.mean():.4f}")
 
 
def quantile_regression_intuition():
    """
    Explain why quantile regression matters.
    """
    print("
 
Quantile Regression: Beyond the Mean")
    print("=" * 60)
    
    print("
Why predict quantiles instead of mean?")
    print("-" * 40)
    
    print("""
    Scenario: Predicting delivery times
    
    • Mean prediction: 5 days
    • Median prediction (Q(0.5)): 4 days  
    • 90th percentile (Q(0.9)): 12 days
    
    A customer asking "when will it arrive?" needs different info:
    - "Most likely": Use median (robust to outliers)
    - "Plan for worst case": Use high quantile (Q(0.9))
    - "Optimistic estimate": Use low quantile (Q(0.1))
    
    Mean-only models can't provide this richness!
    """)
    
    # Simulate heteroscedastic data
    np.random.seed(42)
    x = np.linspace(0, 10, 1000)
    # Variance increases with x (heteroscedasticity)
    y = 2*x + 1 + np.random.normal(0, 1 + 0.5*x, 1000)
    
    # Quantiles at different x values
    print("Simulated delivery time quantiles:")
    print(f"{'x (distance)':<15} {'Q(0.1)':<10} {'Q(0.5)':<10} {'Q(0.9)':<10}")
    print("-" * 45)
    
    for x_val in [1, 5, 10]:
        mask = np.abs(x - x_val) < 0.5
        local_y = y[mask]
        q10 = np.percentile(local_y, 10)
        q50 = np.percentile(local_y, 50)
        q90 = np.percentile(local_y, 90)
        print(f"{x_val:<15.0f} {q10:<10.2f} {q50:<10.2f} {q90:<10.2f}")
    
    print("
Note: Spread (Q90 - Q10) increases with distance!")
 
 
demonstrate_quantiles()
inverse_transform_sampling()
quantile_regression_intuition()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134

import numpy as np
from scipy import stats
 
def ks_test_demonstration():
    """
    Demonstrate the Kolmogorov-Smirnov test for distribution comparison.
    """
    print("Kolmogorov-Smirnov Test")
    print("=" * 60)
    
    np.random.seed(42)
    
    # Test 1: Data from N(0,1) vs N(0,1) hypothesis (should match)
    print("
Test 1: N(0,1) data vs N(0,1) hypothesis")
    data_1 = np.random.normal(0, 1, 500)
    stat_1, pval_1 = stats.kstest(data_1, 'norm', args=(0, 1))
    print(f"  KS statistic: {stat_1:.4f}")
    print(f"  p-value: {pval_1:.4f}")
    print(f"  Conclusion: {'Match (p > 0.05)' if pval_1 > 0.05 else 'Mismatch'}")
    
    # Test 2: Data from N(0.5, 1) vs N(0,1) hypothesis (should differ)
    print("
Test 2: N(0.5, 1) data vs N(0,1) hypothesis")
    data_2 = np.random.normal(0.5, 1, 500)
    stat_2, pval_2 = stats.kstest(data_2, 'norm', args=(0, 1))
    print(f"  KS statistic: {stat_2:.4f}")
    print(f"  p-value: {pval_2:.4f}")
    print(f"  Conclusion: {'Match (p > 0.05)' if pval_2 > 0.05 else 'Mismatch detected!'}")
    
    # Test 3: Two-sample KS test
    print("
Test 3: Two-sample KS test")
    sample_a = np.random.normal(0, 1, 300)
    sample_b = np.random.normal(0, 1.2, 300)  # Slightly different variance
    
    stat_3, pval_3 = stats.ks_2samp(sample_a, sample_b)
    print(f"  Sample A: N(0, 1), n=300")
    print(f"  Sample B: N(0, 1.44), n=300")
    print(f"  KS statistic: {stat_3:.4f}")
    print(f"  p-value: {pval_3:.4f}")
    print(f"  Conclusion: {'Same distribution' if pval_3 > 0.05 else 'Different distributions!'}")
 
 
def empirical_cdf():
    """
    Demonstrate empirical CDF construction.
    """
    print("
 
Empirical CDF (ECDF)")
    print("=" * 60)
    
    print("
The ECDF is the ML estimator for the true CDF:")
    print("  F_n(x) = (1/n) * #{i : X_i ≤ x}")
    print("  = Fraction of samples at or below x")
    
    np.random.seed(42)
    samples = np.random.exponential(scale=2, size=20)
    samples_sorted = np.sort(samples)
    
    print(f"
Samples (sorted): {samples_sorted[:8].round(2)}...")
    
    print("
ECDF evaluation:")
    for x in [0.5, 1.0, 2.0, 5.0]:
        ecdf_val = np.mean(samples <= x)
        true_cdf = stats.expon(scale=2).cdf(x)
        print(f"  F̂({x:.1f}) = {ecdf_val:.3f}, True F({x:.1f}) = {true_cdf:.3f}")
    
    print("
Glivenko-Cantelli Theorem:")
    print("  As n → ∞, sup|F_n(x) - F(x)| → 0 almost surely")
    print("  The ECDF converges uniformly to the true CDF!")
    
    # Show convergence
    print("
Convergence with sample size:")
    for n in [10, 100, 1000, 10000]:
        samples_n = np.random.exponential(scale=2, size=n)
        # KS statistic measures sup|ECDF - CDF|
        ks_stat, _ = stats.kstest(samples_n, 'expon', args=(0, 2))
        print(f"  n = {n:>5}: max|F̂ - F| ≈ {ks_stat:.4f}")
 
 
def calibration_check():
    """
    Show how CDF/quantiles are used for calibration checking.
    """
    print("
 
Model Calibration via CDF")
    print("=" * 60)
    
    print("
A well-calibrated probabilistic model satisfies:")
    print("  For predicted CDF F, observed y should have F(y) ~ Uniform(0,1)")
    print("
This is the 'PIT' (Probability Integral Transform)")
    
    np.random.seed(42)
    n = 1000
    
    # Scenario 1: Well-calibrated model
    # True: Y ~ N(μ, 1), Model predicts N(μ, 1)
    true_mu = np.random.uniform(-1, 1, n)
    y_true = true_mu + np.random.normal(0, 1, n)
    
    # Model predicts N(true_mu, 1) - correctly specified
    pit_values = stats.norm.cdf(y_true, loc=true_mu, scale=1)
    
    # PIT should be uniform
    ks_stat, p_val = stats.kstest(pit_values, 'uniform')
    print(f"
Well-calibrated model:")
    print(f"  PIT ~ Uniform test: KS = {ks_stat:.4f}, p = {p_val:.4f}")
    print(f"  Calibrated: {'Yes' if p_val > 0.05 else 'No'}")
    
    # Scenario 2: Miscalibrated (underestimated variance)
    # Model thinks variance is 0.5, but true variance is 1
    pit_miscal = stats.norm.cdf(y_true, loc=true_mu, scale=0.5)
    
    ks_stat2, p_val2 = stats.kstest(pit_miscal, 'uniform')
    print(f"
Miscalibrated model (underestimated σ):")
    print(f"  PIT ~ Uniform test: KS = {ks_stat2:.4f}, p = {p_val2:.4f}")
    print(f"  Calibrated: {'Yes' if p_val2 > 0.05 else 'No'}")
 
 
ks_test_demonstration()
empirical_cdf()
calibration_check()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121

import numpy as np
from scipy import stats
from scipy.special import logsumexp
 
def log_space_computations():
    """
    Demonstrate log-space probability computations.
    """
    print("Log-Space Probability Computations")
    print("=" * 60)
    
    # Problem: Compute product of many small probabilities
    print("
Problem: Compute ∏ P(xᵢ) for many observations")
    
    # Example: 1000 samples from N(0,1), compute joint probability
    np.random.seed(42)
    samples = np.random.normal(0, 1, 1000)
    normal = stats.norm(0, 1)
    
    # BAD: Direct product (will underflow)
    probs = normal.pdf(samples)
    direct_product = np.prod(probs)
    print(f"
  Direct product: {direct_product}")  # Will be 0.0
    
    # GOOD: Log space computation
    log_probs = normal.logpdf(samples)
    log_product = np.sum(log_probs)
    print(f"  Log of product (sum of logs): {log_product:.4f}")
    
    # The actual product is exp(-1416) ≈ 10^(-615), way smaller than 10^(-308)
    print(f"  This is approximately 10^{log_product/np.log(10):.0f}")
    print(f"  Double precision limit: ~10^(-308)")
 
 
def logsumexp_trick():
    """
    Demonstrate the log-sum-exp trick for stable probability addition.
    """
    print("
 
Log-Sum-Exp Trick")
    print("=" * 60)
    
    print("
Problem: Compute log(p₁ + p₂) given log(p₁), log(p₂)")
    print("Naive: log(exp(log_p1) + exp(log_p2)) → overflow/underflow")
    
    # Example with very small probabilities
    log_p1 = -500  # p1 ≈ 10^(-217)
    log_p2 = -501  # p2 ≈ 10^(-218)
    
    print(f"
  log(p₁) = {log_p1}, log(p₂) = {log_p2}")
    
    # NAive approach fails
    try:
        naive_result = np.log(np.exp(log_p1) + np.exp(log_p2))
        print(f"  Naive result: {naive_result}")
    except:
        print("  Naive: OVERFLOW/UNDERFLOW")
    
    # Actually np.exp(-500) = 0 due to underflow
    print(f"  np.exp({log_p1}) = {np.exp(log_p1)}")  # 0.0
    print(f"  So naive gives: log(0 + 0) = -inf")
    
    # Log-sum-exp trick
    # log(exp(a) + exp(b)) = a + log(1 + exp(b-a))  [if b < a]
    stable_result = logsumexp([log_p1, log_p2])
    print(f"
  Log-sum-exp result: {stable_result:.6f}")
    
    # Verify: log(p1 + p2) when p1/p2 = e^1, so p1 ≈ 2.7 * p2
    # log(p1 + p2) = log(p1 * (1 + p2/p1)) = log_p1 + log(1 + e^(-1))
    expected = log_p1 + np.log(1 + np.exp(log_p2 - log_p1))
    print(f"  Expected (analytical): {expected:.6f}")
 
 
def extreme_tail_probabilities():
    """
    Demonstrate handling of extreme tail probabilities.
    """
    print("
 
Extreme Tail Probabilities")
    print("=" * 60)
    
    normal = stats.norm(0, 1)
    
    print("
P(X > x) for large x in N(0,1):")
    print(f"{'x':>8} {'1-cdf(x)':>20} {'sf(x)':>20} {'logsf(x)':>15}")
    print("-" * 65)
    
    for x in [3, 5, 8, 10, 30, 37]:
        one_minus_cdf = 1 - normal.cdf(x)
        survival = normal.sf(x)
        log_survival = normal.logsf(x)
        
        print(f"{x:>8} {one_minus_cdf:>20.10e} {survival:>20.10e} {log_survival:>15.4f}")
    
    print("
Note: 1 - cdf(x) loses precision for large x")
    print("      sf(x) and logsf(x) remain accurate")
    
    # Critical for p-value computations
    print("
 
This matters for p-values:")
    z_score = 8.0
    print(f"  z-score = {z_score}")
    print(f"  Two-tailed p = 2 * P(|X| > {z_score})")
    print(f"  Using 1 - cdf: {2 * (1 - normal.cdf(z_score)):.2e}")
    print(f"  Using sf:      {2 * normal.sf(z_score):.2e}")
    print(f"  Using logsf:   exp(log(2) + logsf) = {np.exp(np.log(2) + normal.logsf(z_score)):.2e}")
 
 
log_space_computations()
logsumexp_trick()
extreme_tail_probabilities()