第7章 信息论与统计学（Information Theory and Statistics）
一、导读
7.1 本章学习目标
通过本章的学习，你将能够：
理解信息论与统计学的深刻联系：掌握KL散度在假设检验中的作用
掌握充分统计量的信息论刻画：理解为什么充分统计量不损失信息
理解Fisher信息与KL散度的关系：掌握Cramér-Rao下界的信息论解释
掌握最大熵原理及其应用：能够在给定约束下推导最大熵分布
理解EM算法的信息论视角：从KL散度最小化角度理解期望最大化
掌握AIC/BIC准则的信息论基础：理解模型选择中的复杂度惩罚
能够应用互信息进行特征选择：理解基于信息论的特征选择方法
理解信息投影的概念：掌握I-投影和M-投影的区别
7.2 为什么本章如此重要
信息论与统计学有着天然的联系。两者都关注概率分布、不确定性、以及从数据中学习。信息论提供了度量分布差异、量化信息含量、评估模型复杂度的工具，这些工具在现代统计学中无处不在。
理论重要性：
KL散度是统计推断的核心工具。它不仅是假设检验中错误概率的指数衰减速率（Chernoff-Stein引理），也是最大似然估计的渐近性质分析的基础。Fisher信息矩阵可以从KL散度的二阶导数得到，这建立了局部几何与信息论的联系。
实践重要性：
在机器学习中，交叉熵损失本质上是KL散度；EM算法可以看作交替最小化KL散度；变分推断的目标函数ELBO包含KL散度项。理解这些信息论概念有助于深入掌握现代机器学习算法。
方法论重要性：
最大熵原理提供了一种系统的建模方法：在满足已知约束的前提下，选择最"无偏"的分布。这种方法在自然语言处理、图像处理、金融建模等领域都有应用。
7.3 核心问题引导
在学习本章之前，请思考以下问题：
KL散度如何度量两个分布的"距离"？ 为什么它不是真正的距离？在假设检验中，KL散度有什么统计意义？
什么是充分统计量？ 为什么说充分统计量包含了关于参数的全部信息？如何用数据处理不等式刻画充分性？
Fisher信息是什么？ 它与参数估计的精度有什么关系？为什么Fisher信息可以从KL散度导出？
最大熵原理的直觉是什么？ 为什么在只知均值和方差的情况下，高斯分布是最合理的假设？
EM算法为什么要引入隐变量？ 从KL散度角度看，E步和M步分别在做什么？
模型选择中为什么要惩罚复杂度？ AIC和BIC的惩罚项有什么信息论解释？
互信息如何用于特征选择？ 与相关系数相比，互信息有什么优势？
什么是I-投影？ 它在指数族分布的理论中起什么作用？
二、核心概念详解
7.1 KL散度与假设检验
7.1.1 二元假设检验
考虑二元假设检验问题：
- $H_0$: 数据来自分布 $P$
- $H_1$: 数据来自分布 $Q$
给定观测 $x^n = (x_1, \ldots, x_n)$，决策规则 $\phi: \mathcal{X}^n \to \{0, 1\}$。
两类错误：
- 第一类错误（虚警）：$\alpha = P(\phi(x^n) = 1)$
- 第二类错误（漏检）：$\beta = Q(\phi(x^n) = 0)$
7.1.2 Chernoff-Stein引理
定理（Chernoff-Stein引理）：在固定第一类错误 $\alpha \leq \epsilon < 1$ 的条件下，第二类错误 $\beta$ 的最优指数衰减速率为：
$$\lim_{n \to \infty} -\frac{1}{n} \log_2 \beta_n^* = D(P || Q)$$
其中 $\beta_n^*$ 是最小可能的第二类错误。
直观理解：KL散度 $D(P||Q)$ 度量了两个分布的"可区分度"。KL散度越大，越容易从样本中区分出两个分布，第二类错误的衰减速率越快。
证明思路：
使用典型集论证。定义 $P$ 的典型集 $A_\epsilon^{(n)}(P)$。决策规则：如果 $x^n \in A_\epsilon^{(n)}(P)$，接受 $H_0$；否则接受 $H_1$。
- 第一类错误：$P(x^n \notin A_\epsilon^{(n)}(P)) < \epsilon$（由AEP）
- 第二类错误：$\beta = Q(x^n \in A_\epsilon^{(n)}(P))$
对于典型集中的序列，$p(x^n) \approx 2^{-nH(P)}$，因此：
$$\beta \approx |A_\epsilon^{(n)}(P)| \cdot 2^{-nH(Q)} \approx 2^{nH(P)} \cdot 2^{-nH(Q)} = 2^{-nD(P||Q)}$$
更严格的证明需要使用类型法（Method of Types）。
代码示例7-1：假设检验错误概率仿真
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
def simulate_hypothesis_test(n, p_dist, q_dist, n_trials=1000):
"""
模拟二元假设检验
参数:
n: 样本数
p_dist: H0下的分布
q_dist: H1下的分布
n_trials: 试验次数
返回:
alpha, beta (两类错误率)
"""
type1_errors = 0
type2_errors = 0
# 阈值：基于似然比
# 简化：使用样本均值作为统计量
for _ in range(n_trials):
# H0为真
samples_p = p_dist.rvs(n)
mean_p = np.mean(samples_p)
# H1为真
samples_q = q_dist.rvs(n)
mean_q = np.mean(samples_q)
# 决策规则：如果样本均值 > 阈值，选择H1
threshold = (p_dist.mean() + q_dist.mean()) / 2
if mean_p > threshold:
type1_errors += 1
if mean_q <= threshold:
type2_errors += 1
alpha = type1_errors / n_trials
beta = type2_errors / n_trials
return alpha, beta
def kl_divergence_normal(p_mean, p_var, q_mean, q_var):
"""高斯分布间的KL散度"""
return 0.5 * (np.log(q_var/p_var) + (p_var + (p_mean-q_mean)**2)/q_var - 1) / np.log(2)
# 测试
p = norm(loc=0, scale=1)
q = norm(loc=1, scale=1)
D_PQ = kl_divergence_normal(0, 1, 1, 1)
print(f"D(P||Q) = {D_PQ:.4f} bits")
n_values = [10, 20, 50, 100, 200]
beta_values = []
for n in n_values:
alpha, beta = simulate_hypothesis_test(n, p, q, n_trials=500)
beta_values.append(beta)
print(f"n={n:3d}: alpha={alpha:.4f}, beta={beta:.4f}, -log(beta)/n={-np.log2(beta+1e-30)/n:.4f}")
# 可视化
plt.figure(figsize=(10, 6))
plt.semilogy(n_values, beta_values, 'bo-', linewidth=2, markersize=8)
plt.axhline(y=2**(-D_PQ), color='r', linestyle='--', alpha=0.5, label=f'2^(-D(P||Q))={2**(-D_PQ):.4f}')
plt.xlabel('样本数 n', fontsize=12)
plt.ylabel('第二类错误 β', fontsize=12)
plt.title('假设检验错误概率随样本数的变化', fontsize=14)
plt.grid(True, alpha=0.3)
plt.legend()
plt.tight_layout()
plt.savefig('hypothesis_test.png', dpi=150)
plt.show()7.1.3 Sanov定理
Sanov定理是Chernoff-Stein引理的推广，处理多个假设的情况。
定理（Sanov定理）：设 $X_1, \ldots, X_n$ i.i.d.服从 $P$。对于分布集合 $E$：
$$P(\hat{P}_n \in E) \approx 2^{-n \min_{Q \in E} D(Q || P)}$$
其中 $\hat{P}_n$ 是经验分布。
应用：大偏差理论、鲁棒统计、信息几何。
7.2 充分统计量
7.2.1 充分性的定义
定义：统计量 $T(X)$ 称为关于参数 $\theta$ 充分的，如果给定 $T(X)$，样本 $X$ 的条件分布不依赖于 $\theta$。
等价地（因子分解定理）：$T(X)$ 是充分的当且仅当似然函数可以分解为：
$$p(x|\theta) = g(T(x), \theta) h(x)$$
7.2.2 信息论刻画
定理：$T(X)$ 是充分统计量当且仅当：
$$I(\theta; X) = I(\theta; T(X))$$
即 $T(X)$ 保留了关于 $\theta$ 的全部信息。
证明：
考虑马尔可夫链 $\theta \to X \to T(X)$。
由数据处理不等式：$I(\theta; X) \geq I(\theta; T(X))$。
等号成立当且仅当 $\theta \to T(X) \to X$ 也形成马尔可夫链，即给定 $T(X)$，$X$ 与 $\theta$ 独立。这正是充分性的定义。
例子：
对于i.i.d.高斯样本 $X_1, \ldots, X_n \sim \mathcal{N}(\mu, 1)$：
- 样本均值 $\bar{X} = \frac{1}{n}\sum X_i$ 是关于 $\mu$ 的充分统计量
- $I(\mu; X^n) = I(\mu; \bar{X})$
7.2.3 最小充分统计量
定义：如果 $T(X)$ 是充分的，且对于任何其他充分统计量 $S(X)$，存在函数 $f$ 使得 $T(X) = f(S(X))$，则 $T(X)$ 是最小充分统计量。
最小充分统计量是数据的最大压缩表示，不损失关于参数的信息。
7.3 Fisher信息与Cramér-Rao下界
7.3.1 Fisher信息的定义
定义：对于参数化分布 $p(x|\theta)$，Fisher信息定义为：
$$I(\theta) = \mathbb{E}\left[\left(\frac{\partial}{\partial \theta} \log p(X|\theta)\right)^2\right]$$
在多参数情况下，Fisher信息矩阵：
$$[I(\theta)]_{ij} = \mathbb{E}\left[\frac{\partial}{\partial \theta_i} \log p(X|\theta) \cdot \frac{\partial}{\partial \theta_j} \log p(X|\theta)\right]$$
7.3.2 与KL散度的关系
定理：KL散度的二阶导数等于Fisher信息：
$$\frac{\partial^2}{\partial \theta'^2} D(p(\cdot|\theta) || p(\cdot|\theta'))\bigg|_{\theta'=\theta} = I(\theta)$$
证明：
$$\begin{align*}
D(p_\theta || p_{\theta'}) &= \int p(x|\theta) \log \frac{p(x|\theta)}{p(x|\theta')} dx \\
\frac{\partial}{\partial \theta'} D &= -\int p(x|\theta) \frac{\partial}{\partial \theta'} \log p(x|\theta') dx \\
\frac{\partial^2}{\partial \theta'^2} D &= -\int p(x|\theta) \frac{\partial^2}{\partial \theta'^2} \log p(x|\theta') dx
\end{align*}$$
在 $\theta' = \theta$ 处：
$$\frac{\partial^2}{\partial \theta'^2} D\bigg|_{\theta'=\theta} = -\int p(x|\theta) \frac{\partial^2}{\partial \theta^2} \log p(x|\theta) dx = I(\theta)$$
（最后一步使用了得分函数的性质）
直观理解：Fisher信息度量了分布对参数变化的敏感度。KL散度的曲率越大，分布随参数变化越快，参数越容易被估计。
7.3.3 Cramér-Rao下界
定理（Cramér-Rao下界）：对于无偏估计器 $\hat{\theta}$：
$$\text{Var}(\hat{\theta}) \geq \frac{1}{I(\theta)}$$
信息论解释：Fisher信息越大，参数估计的最小方差越小。这与KL散度的解释一致：分布对参数敏感 → 容易估计 → 方差小。
代码示例7-2：Fisher信息计算
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
def fisher_info_normal_mean(n, sigma):
"""
高斯分布均值的Fisher信息
对于 N(μ, σ²)，关于μ的Fisher信息为 n/σ²
"""
return n / sigma**2
def cramer_rao_bound(n, sigma):
"""Cramér-Rao下界"""
return sigma**2 / n
# 测试
sigma = 1.0
n_values = np.arange(1, 101)
fi_values = [fisher_info_normal_mean(n, sigma) for n in n_values]
crb_values = [cramer_rao_bound(n, sigma) for n in n_values]
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))
ax1.plot(n_values, fi_values, 'b-', linewidth=2)
ax1.set_xlabel('样本数 n', fontsize=12)
ax1.set_ylabel('Fisher信息', fontsize=12)
ax1.set_title(f'Fisher信息 vs 样本数 (σ={sigma})', fontsize=14)
ax1.grid(True, alpha=0.3)
ax2.plot(n_values, crb_values, 'r-', linewidth=2)
ax2.set_xlabel('样本数 n', fontsize=12)
ax2.set_ylabel('方差下界', fontsize=12)
ax2.set_title(f'Cramér-Rao下界 vs 样本数 (σ={sigma})', fontsize=14)
ax2.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('fisher_info.png', dpi=150)
plt.show()
print("Fisher信息与Cramér-Rao下界")
print("=" * 60)
for n in [1, 10, 100]:
fi = fisher_info_normal_mean(n, sigma)
crb = cramer_rao_bound(n, sigma)
print(f"n={n:3d}: Fisher信息={fi:.4f}, CR下界={crb:.4f}")7.4 最大熵原理
7.4.1 原理陈述
最大熵原理：在给定约束条件下，选择具有最大熵的概率分布。
直觉：最大熵分布是对已知信息的最"保守"利用——它满足所有已知约束，但对未知部分不做任何额外假设。
7.4.2 常见约束下的最大熵分布
| 约束 | 最大熵分布 |
|---|---|
| 无约束（有限字母表） | 均匀分布 |
| 均值约束 | 指数分布 |
| 均值和方差约束 | 高斯分布 |
| 正半轴 + 均值约束 | 指数分布 |
| 离散 + 均值约束 | 几何分布 |
7.4.3 推导方法
使用拉格朗日乘子法。以均值和方差约束为例：
最大化 $H(p) = -\int p(x) \log p(x) dx$
约束：
$\int p(x) dx = 1$
$\int x p(x) dx = \mu$
$\int (x-\mu)^2 p(x) dx = \sigma^2$
拉格朗日函数：
$$L = -\int p \log p + \lambda_0(\int p - 1) + \lambda_1(\int xp - \mu) + \lambda_2(\int (x-\mu)^2 p - \sigma^2)$$
求变分导数并令为零：
$$\frac{\delta L}{\delta p} = -\log p - 1 + \lambda_0 + \lambda_1 x + \lambda_2 (x-\mu)^2 = 0$$
解得：
$$p(x) \propto \exp(\lambda_1 x + \lambda_2 (x-\mu)^2)$$
代入约束条件，可得高斯分布。
代码示例7-3：最大熵分布验证
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm, expon, uniform
def verify_max_entropy():
"""验证不同约束下的最大熵分布"""
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
# 1. 无约束（有限区间）→ 均匀分布
x1 = np.linspace(0, 1, 100)
uniform_dist = uniform.pdf(x1, 0, 1)
axes[0, 0].plot(x1, uniform_dist, 'b-', linewidth=2)
axes[0, 0].set_title('无约束 → 均匀分布', fontsize=12)
axes[0, 0].set_ylim(0, 1.5)
axes[0, 0].grid(True, alpha=0.3)
# 2. 正半轴 + 均值 → 指数分布
x2 = np.linspace(0, 5, 100)
lambda_param = 1.0 # 均值 = 1/lambda = 1
exp_dist = expon.pdf(x2, scale=1/lambda_param)
axes[0, 1].plot(x2, exp_dist, 'r-', linewidth=2)
axes[0, 1].set_title('正半轴 + 均值约束 → 指数分布', fontsize=12)
axes[0, 1].grid(True, alpha=0.3)
# 3. 均值 + 方差 → 高斯分布
x3 = np.linspace(-4, 4, 100)
mu, sigma = 0, 1
gauss_dist = norm.pdf(x3, mu, sigma)
axes[1, 0].plot(x3, gauss_dist, 'g-', linewidth=2)
axes[1, 0].set_title('均值 + 方差约束 → 高斯分布', fontsize=12)
axes[1, 0].grid(True, alpha=0.3)
# 4. 比较熵值
distributions = ['均匀', '指数', '高斯']
# 均匀分布熵
a, b = 0, 1
H_uniform = np.log(b - a)
# 指数分布熵
H_exp = 1 - np.log(lambda_param)
# 高斯分布熵
H_gauss = 0.5 * np.log(2 * np.pi * np.e * sigma**2)
entropies = [H_uniform, H_exp, H_gauss]
axes[1, 1].bar(distributions, entropies, color=['blue', 'red', 'green'], alpha=0.7)
axes[1, 1].set_ylabel('微分熵 (nats)', fontsize=12)
axes[1, 1].set_title('各最大熵分布的熵值', fontsize=12)
axes[1, 1].grid(True, alpha=0.3, axis='y')
for i, v in enumerate(entropies):
axes[1, 1].text(i, v + 0.1, f'{v:.4f}', ha='center')
plt.tight_layout()
plt.savefig('max_entropy_distributions.png', dpi=150)
plt.show()
print("最大熵分布及其熵值")
print("=" * 40)
print(f"均匀分布 U(0,1): H = {H_uniform:.4f} nats")
print(f"指数分布 Exp(1): H = {H_exp:.4f} nats")
print(f"高斯分布 N(0,1): H = {H_gauss:.4f} nats")
verify_max_entropy()7.4.4 应用
自然语言处理：最大熵模型用于词性标注、命名实体识别等任务。
生态学：最大熵方法用于物种分布建模。
金融：在已知某些矩约束下，用最大熵分布建模资产收益。
7.5 EM算法的信息论视角
7.5.1 EM算法回顾
EM算法用于含有隐变量的概率模型参数估计。
E步：计算隐变量的后验期望
$$Q(\theta|\theta^{(t)}) = \mathbb{E}_{Z|X,\theta^{(t)}}[\log p(X,Z|\theta)]$$
M步：最大化期望
$$\theta^{(t+1)} = \arg\max_\theta Q(\theta|\theta^{(t)})$$
7.5.2 KL散度解释
EM算法可以看作交替最小化以下KL散度：
定义辅助分布 $q(z)$，目标是最小化：
$$D(q(z) || p(z|x,\theta))$$
E步：固定 $\theta$，最优 $q(z) = p(z|x,\theta^{(t)})$
M步：固定 $q$，最大化关于 $\theta$ 的期望对数似然
这等价于最大化ELBO（Evidence Lower Bound）：
$$\text{ELBO}(q, \theta) = \mathbb{E}_q[\log p(X,Z|\theta)] - \mathbb{E}_q[\log q(Z)]$$
因为：
$$\log p(X|\theta) = \text{ELBO}(q, \theta) + D(q(Z) || p(Z|X,\theta)) \geq \text{ELBO}(q, \theta)$$
7.5.3 收敛性
EM算法保证ELBO单调增加，因此对数似然也单调增加（或至少不减少）。
7.6 模型选择准则
7.6.1 AIC（Akaike Information Criterion）
$$\text{AIC} = -2 \log L(\hat{\theta}) + 2k$$
其中 $k$ 是参数个数。
信息论解释：AIC估计了模型与真实分布之间的KL散度。第一项衡量拟合优度，第二项惩罚复杂度。
7.6.2 BIC（Bayesian Information Criterion）
$$\text{BIC} = -2 \log L(\hat{\theta}) + k \log n$$
BIC的惩罚项更强，倾向于选择更简单的模型。
推导：BIC近似于贝叶斯模型证据（marginal likelihood）。
7.6.3 MDL（Minimum Description Length）
MDL原则选择使数据和模型的总描述长度最短的模型。
与信息论的联系：描述长度对应于编码成本，最小化描述长度等价于最大化压缩率。
代码示例7-4：模型选择比较
import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures
def compute_aic(n, k, rss):
"""计算AIC"""
sigma2_hat = rss / n
return n * np.log(sigma2_hat) + 2 * k
def compute_bic(n, k, rss):
"""计算BIC"""
sigma2_hat = rss / n
return n * np.log(sigma2_hat) + k * np.log(n)
# 生成数据
np.random.seed(42)
n = 100
x = np.linspace(0, 10, n).reshape(-1, 1)
y_true = 2 * x + 0.5 * x**2
noise = np.random.normal(0, 2, size=n).reshape(-1, 1)
y = y_true + noise
# 尝试不同阶数的多项式
degrees = range(1, 11)
aic_values = []
bic_values = []
rss_values = []
for d in degrees:
poly = PolynomialFeatures(degree=d)
x_poly = poly.fit_transform(x)
model = LinearRegression()
model.fit(x_poly, y)
y_pred = model.predict(x_poly)
rss = np.sum((y - y_pred)**2)
k = d + 1 # 参数个数（包括截距）
aic = compute_aic(n, k, rss)
bic = compute_bic(n, k, rss)
aic_values.append(aic)
bic_values.append(bic)
rss_values.append(rss)
# 可视化
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))
ax1.plot(degrees, rss_values, 'bo-', linewidth=2, markersize=8, label='RSS')
ax1.set_xlabel('多项式阶数', fontsize=12)
ax1.set_ylabel('残差平方和', fontsize=12)
ax1.set_title('拟合优度 vs 模型复杂度', fontsize=14)
ax1.grid(True, alpha=0.3)
ax1.legend()
ax2.plot(degrees, aic_values, 'rs-', linewidth=2, markersize=8, label='AIC')
ax2.plot(degrees, bic_values, 'g^-', linewidth=2, markersize=8, label='BIC')
ax2.set_xlabel('多项式阶数', fontsize=12)
ax2.set_ylabel('准则值', fontsize=12)
ax2.set_title('AIC/BIC vs 模型复杂度', fontsize=14)
ax2.grid(True, alpha=0.3)
ax2.legend()
plt.tight_layout()
plt.savefig('model_selection.png', dpi=150)
plt.show()
print("模型选择结果")
print("=" * 60)
best_aic_deg = degrees[np.argmin(aic_values)]
best_bic_deg = degrees[np.argmin(bic_values)]
print(f"AIC最佳阶数: {best_aic_deg}")
print(f"BIC最佳阶数: {best_bic_deg}")
print(f"真实阶数: 2")7.7 互信息与特征选择
7.7.1 基于互信息的特征选择
思想：选择与标签互信息最大的特征。
$$\text{Score}(X_i) = I(X_i; Y)$$
优势：
- 捕捉非线性关系（相关系数只能检测线性关系）
- 适用于混合数据类型（连续+离散）
挑战：
- 估计互信息需要足够的样本
- 高维情况下估计困难
7.7.2 mRMR（最小冗余最大相关）
同时考虑特征与标签的相关性和特征之间的冗余：
$$\text{Score}(X_i) = I(X_i; Y) - \frac{1}{|S|}\sum_{X_j \in S} I(X_i; X_j)$$
其中 $S$ 是已选特征集合。
代码示例7-5：互信息特征选择
import numpy as np
from sklearn.datasets import load_iris
from sklearn.feature_selection import mutual_info_classif
# 加载数据
iris = load_iris()
X, y = iris.data, iris.target
# 计算互信息
mi_scores = mutual_info_classif(X, y, discrete_features=False, random_state=42)
print("Iris数据集特征的互信息得分")
print("=" * 50)
for name, score in zip(iris.feature_names, mi_scores):
print(f"{name:25s}: {score:.4f}")
# 排序
sorted_indices = np.argsort(mi_scores)[::-1]
print("\n特征重要性排名:")
for rank, idx in enumerate(sorted_indices, 1):
print(f" {rank}. {iris.feature_names[idx]}: {mi_scores[idx]:.4f}")三、程序员的视角
7.1 在机器学习中的应用
交叉熵损失：分类任务中最常用的损失函数，本质上是KL散度加上常数。
变分推断：通过最小化KL散度来近似后验分布。
GAN：生成对抗网络可以看作最小化JS散度（一种对称化的KL散度）。
7.2 实用建议
估计互信息：使用KSG估计器或基于k近邻的方法
模型选择：小样本用BIC，大样本用AIC
最大熵建模：使用现成库（如MaxEnt）
7.3 常见陷阱
KL散度不对称：$D(p||q) \neq D(q||p)$，选择方向很重要
零概率问题：当 $q(x)=0$ 但 $p(x)>0$ 时，KL散度无穷大
高维估计偏差：互信息估计在高维下有较大偏差
四、代码实战
7.1 完整的信息论统计分析工具箱
import numpy as np
from scipy.stats import entropy as scipy_entropy
from scipy.special import digamma
class InformationStatistics:
"""信息论统计分析工具箱"""
@staticmethod
def kl_divergence_empirical(p_samples, q_pdf):
"""
从样本估计KL散度
参数:
p_samples: 来自分布P的样本
q_pdf: 分布Q的概率密度函数
返回:
D(P||Q)估计
"""
n = len(p_samples)
log_ratios = np.log(q_pdf(p_samples) + 1e-30)
# 使用核密度估计P
from scipy.stats import gaussian_kde
kde = gaussian_kde(p_samples)
log_p = np.log(kde(p_samples) + 1e-30)
return np.mean(log_p - log_ratios)
@staticmethod
def mutual_information_knn(X, Y, k=5):
"""
使用k近邻估计互信息
参数:
X, Y: 数据数组
k: 近邻数
返回:
I(X;Y)估计
"""
from sklearn.neighbors import NearestNeighbors
n = len(X)
# 联合空间
XY = np.column_stack([X, Y])
# 计算距离
nn_xy = NearestNeighbors(n_neighbors=k+1).fit(XY)
dist_xy, _ = nn_xy.kneighbors(XY)
nn_x = NearestNeighbors(n_neighbors=k+1).fit(X.reshape(-1, 1))
dist_x, _ = nn_x.kneighbors(X.reshape(-1, 1))
nn_y = NearestNeighbors(n_neighbors=k+1).fit(Y.reshape(-1, 1))
dist_y, _ = nn_y.kneighbors(Y.reshape(-1, 1))
# KSG估计器
psi_k = digamma(k)
psi_n = digamma(n)
mi = psi_k - psi_n + np.mean(
np.log(dist_x[:, -1] * dist_y[:, -1] / dist_xy[:, -1]**2)
) * (X.shape[1] + Y.shape[1] - X.shape[1] - Y.shape[1]) / 2
return max(0, mi) # 确保非负
@staticmethod
def fisher_info_sample(samples, param_func):
"""
从样本估计Fisher信息
参数:
samples: 样本
param_func: 参数化函数
返回:
Fisher信息估计
"""
n = len(samples)
theta_hat = param_func(samples)
# 数值计算得分函数
epsilon = 1e-5
scores = []
for i in range(n):
# 留一法
samples_minus = np.delete(samples, i)
theta_minus = param_func(samples_minus)
scores.append((theta_hat - theta_minus) / (1/n))
return np.var(scores)
# 使用示例
if __name__ == "__main__":
stats = InformationStatistics()
# KL散度估计
p_samples = np.random.normal(0, 1, 1000)
from scipy.stats import norm
D = stats.kl_divergence_empirical(p_samples, lambda x: norm.pdf(x, 0, 1))
print(f"D(P||Q)估计: {D:.4f}")
# 互信息估计
X = np.random.randn(500)
Y = X + np.random.randn(500) * 0.5
mi = stats.mutual_information_knn(X, Y)
print(f"I(X;Y)估计: {mi:.4f}")五、深入思考
5.1 历史背景
最大熵原理由Jaynes在1957年提出，将信息论与统计力学联系起来。AIC由Akaike在1974年提出，开创了信息准则的先河。EM算法由Dempster, Laird, Rubin在1977年系统阐述，但类似思想早已出现。
5.2 前沿方向
信息几何：将统计流形视为黎曼流形，Fisher信息作为度量张量。
率失真感知理论：结合率失真理论和感知质量度量。
信息瓶颈：在深度学习中平衡压缩和预测。
5.3 延伸阅读
Cover & Thomas《信息论基础》第11章
Jaynes "Information Theory and Statistical Mechanics"
Bishop《Pattern Recognition and Machine Learning》第9章（EM算法）
六、实验笔记
实验1：最大熵分布拟合
目标：验证给定约束下，最大熵分布确实是最优选择。
步骤：
生成具有特定均值和方差的数据
拟合高斯分布、均匀分布、指数分布等
计算各分布与经验分布的KL散度
验证高斯分布的KL散度最小
实验2：互信息特征选择
目标：比较基于互信息和基于相关系数的特征选择方法。
步骤：
生成包含线性和非线性关系的数据
分别使用互信息和相关系数评分特征
训练分类器，比较两种方法选择的特征子集的性能
七、本章小结
本章探讨了信息论与统计学的深刻联系。
KL散度是假设检验中错误概率衰减速率的度量，也是模型选择和信息几何的核心工具。
充分统计量保留了关于参数的全部信息，可以用数据处理不等式精确刻画。
Fisher信息从KL散度的二阶导数导出，决定了参数估计的最小方差（Cramér-Rao下界）。
最大熵原理提供了一种系统的建模方法，在给定约束下选择最无偏的分布。
EM算法可以从KL散度最小化的角度理解，E步和M步交替优化辅助分布和参数。
模型选择准则（AIC/BIC）有明确的信息论解释，平衡拟合优度和模型复杂度。
互信息是强大的特征选择工具，能够捕捉非线性关系。
八、习题
习题7.1（KL散度）：计算两个高斯分布 $\mathcal{N}(\mu_1, \sigma_1^2)$ 和 $\mathcal{N}(\mu_2, \sigma_2^2)$ 之间的KL散度。
习题7.2（充分统计量）：证明对于泊松分布 $Poisson(\lambda)$，样本和 $\sum X_i$ 是充分统计量。
习题7.3（Fisher信息）：计算伯努利分布 $Bernoulli(p)$ 的Fisher信息。
习题7.4（最大熵）：推导在正半轴上、给定均值约束下的最大熵分布。
习题7.5（编程题）：实现基于互信息的特征选择算法，在UCI数据集上测试其效果。
习题7.6（EM算法）：使用EM算法拟合高斯混合模型，可视化每次迭代的ELBO变化。
习题7.7（AIC/BIC）：在线性回归中，比较AIC、BIC和交叉验证选择的模型阶数。
习题7.8（挑战题）：证明Sanov定理的特殊情况（二元假设检验）。
习题7.9（应用题）：在自然语言处理中，如何使用最大熵模型进行文本分类？
习题7.10（开放性问题）：研究信息瓶颈理论在深度学习中的应用，撰写一篇500字左右的分析报告。