Fisher信息量实战:用Python验证Cramér-Rao下界与MLE效率
Fisher信息量实战:用Python验证Cramér-Rao下界与MLE效率
1. 理解Fisher信息量的统计意义
Fisher信息量是衡量概率分布中参数所携带信息多少的重要指标。想象你是一位侦探,Fisher信息量就是你手中的"信息放大镜"——它决定了你从数据中推断参数的精确程度。
对于一个参数θ,Fisher信息量I(θ)定义为:
import numpy as np def fisher_information(p): """计算Bernoulli分布的Fisher信息量""" return 1 / (p * (1 - p))这个简单的公式背后蕴含着深刻的统计意义:
- 当p接近0或1时,Fisher信息量急剧增大,意味着我们更容易准确估计极端概率
- 在p=0.5时信息量最小,此时区分正负类最具挑战性
Fisher信息量的核心性质:
- 与得分函数(对数似然的梯度)的方差直接相关
- 决定了Cramér-Rao下界——任何无偏估计量方差的理论下限
- 与MLE的渐进方差成反比关系
注意:Fisher信息量越大,意味着参数估计可以达到的精度越高,这就像拥有更高像素的相机能捕捉更清晰的图像。
2. Cramér-Rao下界的数学本质
Cramér-Rao不等式为无偏估计量的方差设定了理论下限:
Var(θ̂) ≥ 1/(nI(θ))这个不等式告诉我们,即使是最优的无偏估计量,其方差也不可能突破这个界限。让我们通过Bernoulli分布的例子来验证这一点。
Bernoulli分布的Cramér-Rao下界计算:
| 参数p | Fisher信息量I(p) | n=100时的下界 |
|---|---|---|
| 0.1 | 11.11 | 0.0009 |
| 0.3 | 4.76 | 0.0021 |
| 0.5 | 4.00 | 0.0025 |
3. 构建蒙特卡洛实验框架
我们将设计一个完整的实验来验证MLE的效率是否达到Cramér-Rao下界。以下是实验的关键步骤:
def monte_carlo_experiment(p, n, num_simulations=10000): """执行蒙特卡洛实验验证MLE性质""" fisher_info = fisher_information(p) cr_lower_bound = 1 / (n * fisher_info) mle_variances = [] for _ in range(num_simulations): samples = np.random.binomial(1, p, size=n) p_hat = np.mean(samples) mle_variances.append((p_hat - p)**2) empirical_variance = np.mean(mle_variances) return { 'CR_lower_bound': cr_lower_bound, 'empirical_variance': empirical_variance, 'ratio': empirical_variance / cr_lower_bound }实验设计要点:
- 固定真实参数p和样本量n
- 重复生成随机样本并计算MLE
- 收集MLE的方差并与理论下界比较
4. 结果可视化与分析
让我们用不同样本量进行实验,并可视化结果:
import matplotlib.pyplot as plt sample_sizes = [10, 50, 100, 500, 1000] p = 0.3 results = [monte_carlo_experiment(p, n) for n in sample_sizes] plt.figure(figsize=(10, 6)) plt.plot(sample_sizes, [r['CR_lower_bound'] for r in results], label='Cramér-Rao下界') plt.plot(sample_sizes, [r['empirical_variance'] for r in results], label='MLE经验方差') plt.xlabel('样本量') plt.ylabel('方差') plt.title('MLE方差与Cramér-Rao下界比较(p=0.3)') plt.legend() plt.grid(True) plt.show()典型实验结果会显示:
- 小样本时MLE方差明显高于下界
- 随着样本量增加,MLE方差逐渐接近并最终达到下界
- 验证了MLE的渐进有效性
5. 渐进正态性的实证检验
MLE的另一个重要性质是渐进正态性:
√n(θ̂ - θ) → N(0, 1/I(θ))
我们可以通过QQ图来验证这一性质:
from scipy import stats def check_asymptotic_normality(p, n, num_simulations=5000): """检验MLE的渐进正态性""" mles = [] for _ in range(num_simulations): samples = np.random.binomial(1, p, size=n) mles.append(np.mean(samples)) standardized = (np.array(mles) - p) * np.sqrt(n * fisher_information(p)) stats.probplot(standardized, dist="norm", plot=plt) plt.title(f'QQ图验证渐进正态性(n={n})') plt.show()当样本量足够大时,QQ图上的点将近似落在一条直线上,直观验证了MLE的渐进正态分布性质。
6. 实际应用中的注意事项
虽然理论性质优美,但在实际应用中需要注意:
小样本问题:
- 当np或n(1-p)较小时,正态近似可能不准确
- 考虑使用Wilson区间等改进方法
边界情况处理:
def safe_mle(samples): """处理全0或全1样本的MLE计算""" p_hat = np.mean(samples) n = len(samples) if p_hat == 0: return 1/(2*n) elif p_hat == 1: return 1 - 1/(2*n) return p_hat- 模型误设影响:
- 当真实分布不是Bernoulli时,MLE可能不再有效
- 可通过稳健标准误等方法减轻影响
7. 扩展到其他分布
虽然我们以Bernoulli分布为例,但这种方法可以推广到其他分布。例如,对于泊松分布:
def poisson_fisher_information(lambd): """泊松分布的Fisher信息量""" return 1 / lambd def poisson_mle(samples): """泊松分布的MLE""" return np.mean(samples)不同分布Fisher信息量对比:
| 分布 | 参数 | Fisher信息量公式 |
|---|---|---|
| Bernoulli | p | 1/[p(1-p)] |
| Poisson | λ | 1/λ |
| Normal(μ已知) | σ² | 1/(2σ⁴) |
| Exponential | λ | 1/λ² |
理解这些差异有助于我们在不同建模场景中选择合适的分布和评估估计精度。
