Cumulative Distribution Functions

1import numpy as np
2from scipy import stats
3import matplotlib.pyplot as plt
4
5# ============================================
6# EXAMPLE 1: CDF Evaluation
7# ============================================
8
9# Standard normal CDF
10normal = stats.norm(loc=0, scale=1)
11
12# P(X ≤ 0) = 0.5 (symmetric around 0)
13print(f"P(X ≤ 0) = {normal.cdf(0):.4f}")  # 0.5000
14
15# P(X ≤ 1.96) ≈ 0.975 (the famous 95% CI bound)
16print(f"P(X ≤ 1.96) = {normal.cdf(1.96):.4f}")  # 0.9750
17
18# P(X ≤ -1) ≈ 0.159
19print(f"P(X ≤ -1) = {normal.cdf(-1):.4f}")  # 0.1587
20
21# ============================================
22# EXAMPLE 2: Interval Probabilities
23# ============================================
24
25# P(-1 ≤ X ≤ 1) = F(1) - F(-1)
26prob_interval = normal.cdf(1) - normal.cdf(-1)
27print(f"P(-1 ≤ X ≤ 1) = {prob_interval:.4f}")  # 0.6827 (the 68-95-99.7 rule!)
28
29# P(X > 2) = 1 - F(2)
30prob_greater = 1 - normal.cdf(2)
31print(f"P(X > 2) = {prob_greater:.4f}")  # 0.0228
32
33# ============================================
34# EXAMPLE 3: Inverse CDF (Quantile Function)
35# ============================================
36
37# Q(0.5) = median = 0 for standard normal
38print(f"Median: Q(0.5) = {normal.ppf(0.5):.4f}")  # 0.0000
39
40# Q(0.975) = 1.96 (upper bound for 95% CI)
41print(f"Q(0.975) = {normal.ppf(0.975):.4f}")  # 1.9600
42
43# 95% confidence interval: [Q(0.025), Q(0.975)]
44ci_lower = normal.ppf(0.025)
45ci_upper = normal.ppf(0.975)
46print(f"95% CI: [{ci_lower:.4f}, {ci_upper:.4f}]")  # [-1.96, 1.96]
47
48# ============================================
49# EXAMPLE 4: Inverse Transform Sampling
50# ============================================
51
52np.random.seed(42)
53
54# Generate uniform random numbers
55n_samples = 10000
56u = np.random.uniform(0, 1, size=n_samples)
57
58# Transform to standard normal using inverse CDF
59samples_normal = normal.ppf(u)
60
61# Verify: mean should be ~0, std should be ~1
62print(f"Sample mean: {samples_normal.mean():.4f}")  # ~0
63print(f"Sample std: {samples_normal.std():.4f}")   # ~1
64
65# This is equivalent to np.random.normal(0, 1, n_samples)!
66
67# ============================================
68# EXAMPLE 5: Empirical CDF
69# ============================================
70
71from statsmodels.distributions.empirical_distribution import ECDF
72
73# Generate some data
74data = np.random.normal(0, 1, 1000)
75
76# Build ECDF
77ecdf = ECDF(data)
78
79# Evaluate ECDF at specific points
80print(f"ECDF(0) = {ecdf(0):.4f}")  # ~0.5
81print(f"ECDF(1.96) = {ecdf(1.96):.4f}")  # ~0.975
82
83# ============================================
84# EXAMPLE 6: Kolmogorov-Smirnov Test
85# ============================================
86
87from scipy.stats import kstest
88
89# Test if data comes from standard normal
90ks_stat, p_value = kstest(data, 'norm')
91print(f"KS statistic: {ks_stat:.4f}")
92print(f"P-value: {p_value:.4f}")
93
94# ============================================
95# EXAMPLE 7: Survival Function
96# ============================================
97
98# Exponential survival function
99exp_dist = stats.expon(scale=1/0.5)  # λ = 0.5
100
101# P(X > 2) using survival function
102print(f"P(X > 2) = {exp_dist.sf(2):.4f}")  # 0.3679
103
104# Verify: sf(x) = 1 - cdf(x)
105print(f"1 - F(2) = {1 - exp_dist.cdf(2):.4f}")  # Same!
106
107# ============================================
108# EXAMPLE 8: VAE Reparameterization Trick
109# ============================================
110
111def reparameterize(mu, log_var):
112    """
113    VAE reparameterization trick.
114    Sample z ~ N(mu, exp(log_var)) using inverse transform idea.
115    """
116    std = np.exp(0.5 * log_var)
117    eps = np.random.normal(0, 1, size=mu.shape)  # Standard normal samples
118    return mu + std * eps  # This IS inverse transform sampling!
119
120# Example: encoder outputs mu=2, log_var=0.5
121mu = np.array([2.0])
122log_var = np.array([0.5])
123z = reparameterize(mu, log_var)
124print(f"Sampled z: {z[0]:.4f}")
125
126# ============================================
127# EXAMPLE 9: Probability Calibration Check
128# ============================================
129
130def reliability_diagram(y_true, y_pred_proba, n_bins=10):
131    """Check if predicted probabilities are calibrated."""
132    bin_edges = np.linspace(0, 1, n_bins + 1)
133    bin_means = []
134    bin_true_fractions = []
135
136    for i in range(n_bins):
137        mask = (y_pred_proba >= bin_edges[i]) & (y_pred_proba < bin_edges[i+1])
138        if mask.sum() > 0:
139            bin_means.append(y_pred_proba[mask].mean())
140            bin_true_fractions.append(y_true[mask].mean())
141
142    return np.array(bin_means), np.array(bin_true_fractions)
143
144# ============================================
145# EXAMPLE 10: A/B Testing with CDF
146# ============================================
147
148# Two-sample t-test p-value uses t-distribution CDF
149from scipy.stats import ttest_ind
150
151group_a = np.random.normal(10.0, 2.0, 100)  # Control
152group_b = np.random.normal(10.5, 2.0, 100)  # Treatment
153
154t_stat, p_value = ttest_ind(group_a, group_b)
155print(f"t-statistic: {t_stat:.4f}")
156print(f"p-value: {p_value:.4f}")
157
158# The p-value is computed using CDF of t-distribution!
159# p-value = 2 * (1 - F(|t|)) for two-tailed test
160t_dist = stats.t(df=198)  # degrees of freedom
161manual_pvalue = 2 * (1 - t_dist.cdf(abs(t_stat)))
162print(f"Manual p-value: {manual_pvalue:.4f}")