Python实战：用NumPy模拟大数定律与中心极限定理（附代码）

概率论中的大数定律和中心极限定理是数据分析的基石，但很多初学者往往被抽象的数学公式劝退。今天我们将用Python和NumPy，通过代码实验让这两个重要概念变得直观可见。不需要复杂的数学推导，只需跟着代码一步步操作，你就能亲眼见证这些统计规律如何在数据中显现。

1. 实验环境搭建与基础概念

在开始实验前，我们需要确保环境配置正确。推荐使用Jupyter Notebook进行交互式操作，它能实时显示图表和计算结果。

首先安装必要的库：

pip install numpy matplotlib

大数定律告诉我们，当实验次数足够多时，事件发生的频率会趋近于其理论概率。而中心极限定理则指出，无论原始数据分布如何，样本均值的分布都会趋近于正态分布。这两个定理为统计学中的参数估计和假设检验提供了理论基础。

让我们先定义几个关键术语：

• 样本均值：一组数据的平均值
• 总体均值：理论上的期望值
• 收敛：随着样本量增大，统计量接近某个固定值的过程

2. 模拟大数定律：硬币投掷实验

我们从一个简单的硬币投掷实验开始。理论上，公平硬币出现正面的概率是0.5。让我们用代码模拟这一过程。

import numpy as np
import matplotlib.pyplot as plt

# 设置随机种子保证结果可复现
np.random.seed(42)

# 模拟从1次到10000次投掷
max_trials = 10000
trials = np.arange(1, max_trials+1)
results = np.random.binomial(1, 0.5, max_trials)
cumulative_means = np.cumsum(results) / trials

# 绘制结果
plt.figure(figsize=(10, 6))
plt.plot(trials, cumulative_means, label='实际频率')
plt.axhline(0.5, color='red', linestyle='--', label='理论概率')
plt.xlabel('试验次数')
plt.ylabel('正面出现频率')
plt.title('大数定律演示：硬币投掷实验')
plt.legend()
plt.grid(True)
plt.show()

这段代码会生成一个图表，展示随着试验次数增加，正面出现的频率如何逐渐稳定在0.5附近。你会注意到：

• 在小样本量时（前几百次），频率波动很大
• 随着试验次数增加，波动逐渐减小
• 最终频率收敛到理论概率0.5

常见问题排查：

• 如果图表没有显示，确保已安装matplotlib并调用了plt.show()
• 结果看起来不够平滑？尝试增加max_trials到更大的值
• 想模拟不公平硬币？修改np.random.binomial中的概率参数

3. 中心极限定理的数值验证

中心极限定理更为神奇，它告诉我们无论原始数据分布如何，样本均值的分布都会趋近于正态分布。让我们用三种不同分布来验证这一点。

3.1 均匀分布的例子

def demonstrate_clt(sample_size=30, num_experiments=1000, distribution='uniform'):
    plt.figure(figsize=(12, 5))
    
    # 原始分布
    plt.subplot(1, 2, 1)
    if distribution == 'uniform':
        data = np.random.uniform(0, 1, 10000)
        plt.hist(data, bins=30, density=True, alpha=0.7)
        plt.title('原始均匀分布')
    elif distribution == 'exponential':
        data = np.random.exponential(1, 10000)
        plt.hist(data, bins=30, density=True, alpha=0.7)
        plt.title('原始指数分布')
    else:  # binomial
        data = np.random.binomial(10, 0.3, 10000)
        plt.hist(data, bins=30, density=True, alpha=0.7)
        plt.title('原始二项分布')
    
    # 样本均值分布
    plt.subplot(1, 2, 2)
    sample_means = []
    for _ in range(num_experiments):
        samples = np.random.choice(data, size=sample_size)
        sample_means.append(np.mean(samples))
    
    plt.hist(sample_means, bins=30, density=True, alpha=0.7)
    plt.title(f'样本均值分布 (n={sample_size})')
    plt.tight_layout()
    plt.show()

# 演示均匀分布的情况
demonstrate_clt(distribution='uniform')

3.2 指数分布的例子

# 演示指数分布的情况
demonstrate_clt(distribution='exponential')

3.3 二项分布的例子

# 演示二项分布的情况
demonstrate_clt(distribution='binomial')

这三个例子展示了无论原始数据是均匀分布、指数分布还是二项分布，样本均值的分布都会随着样本量增大而趋近于正态分布。你可以尝试：

1. 修改sample_size参数，观察不同样本量下的分布形态
2. 增加num_experiments使结果更平滑
3. 尝试其他分布类型

4. 进阶应用：置信区间的模拟

理解了中心极限定理后，我们可以用它来构建置信区间。让我们模拟95%置信区间的计算过程。

def simulate_confidence_intervals(population_mean=5, population_std=2, 
                                 sample_size=30, num_simulations=100):
    # 生成正态分布的总体数据
    population = np.random.normal(population_mean, population_std, 10000)
    
    # 存储每次模拟的置信区间
    intervals = []
    contains_mean = []
    
    for _ in range(num_simulations):
        # 从总体中抽样
        sample = np.random.choice(population, size=sample_size)
        sample_mean = np.mean(sample)
        sample_std = np.std(sample, ddof=1)  # 使用无偏估计
        
        # 计算95%置信区间
        margin_of_error = 1.96 * (sample_std / np.sqrt(sample_size))
        ci_lower = sample_mean - margin_of_error
        ci_upper = sample_mean + margin_of_error
        
        intervals.append((ci_lower, ci_upper))
        contains_mean.append(ci_lower <= population_mean <= ci_upper)
    
    # 绘制结果
    plt.figure(figsize=(10, 6))
    for i, (lower, upper) in enumerate(intervals):
        plt.plot([lower, upper], [i, i], color='blue' if contains_mean[i] else 'red')
    plt.axvline(population_mean, color='black', linestyle='--')
    plt.xlabel('数值')
    plt.ylabel('模拟次数')
    plt.title('95%置信区间模拟 (黑色虚线为真实均值)')
    plt.show()
    
    # 计算包含真实均值的比例
    coverage_rate = np.mean(contains_mean)
    print(f"置信区间包含真实均值的比例: {coverage_rate:.2%}")

# 运行模拟
simulate_confidence_intervals()

这个模拟展示了：

• 如何从样本数据构建置信区间
• 为什么95%的置信区间大约有95%的概率包含真实均值
• 当置信区间不包含真实均值时（红色线段），发生了什么

5. 实际应用中的注意事项

虽然这些定理非常强大，但在实际应用中仍需注意以下几点：

1. 样本量要求：

   • 大数定律需要足够多的观测才能显现
   • 中心极限定理通常要求样本量n≥30才能较好近似

2. 数据独立性：

   • 定理假设样本是独立同分布的
   • 现实中很多数据存在自相关（如时间序列数据），这会破坏假设

3. 极端分布：

   • 对于极端偏态或重尾分布，可能需要更大样本量
   • 下表比较了不同分布对中心极限定理收敛速度的影响：

分布类型	收敛速度	建议最小样本量
对称分布（如正态）	快	15-20
中等偏态	中等	30-50
极端偏态	慢	50+

4. 常见误区：
   • 误认为大数定律保证短期结果会"平衡"（赌徒谬误）
   • 忽视中心极限定理对样本量的要求
   • 在明显非独立数据中错误应用这些定理

# 展示样本量不足时的风险
demonstrate_clt(sample_size=5, num_experiments=1000, distribution='exponential')

这个小样本例子清楚地展示了当样本量不足时，样本均值的分布可能还远未接近正态分布。