9.1 动手称量：Python 信息论工具箱

9.1 动手称量：Python 信息论工具箱

本节摘要：交付三件最常用量具的可执行版本——从样本估熵的插件法（附低估偏差的实测数字与校正方向）、互信息的矩阵化用法、AWGN 容量曲线的数值计算；说明工具链的组织方式与"先估后校"的工作习惯。

工具箱的第一课：你算的是估计值

前面各章的代码都在"真值已知"的假设下演算——概率表给定，熵与互信息直接代入。真实工作恰好相反：手上只有样本，分布要靠估。从样本频率直接代入熵公式（插件法）是最直觉的做法，也是偏差的来源：样本有限时，一部分符号恰好没出现、出现的符号频率也被噪声推高或压低，净效果是系统性低估熵——低估的量级随符号数与样本量的比例变化，符号越多、样本越少，偏差越狠。

# 插件法估熵的偏差实测：小样本严重低估
import random
from math import log2

probs = [0.40, 0.20, 0.15, 0.10, 0.06, 0.04, 0.03, 0.02]
H_true = -sum(p * log2(p) for p in probs)

def plugin_entropy(samples, k=8):
    counts = [0] * k
    for s in samples:
        counts[s] += 1
    n = len(samples)
    return -sum((c / n) * log2(c / n) for c in counts if c > 0)

random.seed(11)
print(f"真实熵 = {H_true:.3f} 比特/符号")
for n in [50, 200, 1000, 5000]:
    ests = [plugin_entropy(random.choices(range(8), weights=probs, k=n))
            for _ in range(50)]                       # 50 次重复取平均
    avg = sum(ests) / 50
    print(f"样本量 {n:>5}: 平均估计 {avg:.3f}  偏低 {H_true - avg:.3f}")
# 输出:
# 真实熵 = 2.430 比特/符号
# 样本量    50: 平均估计 2.317  偏低 0.113
# 样本量   200: 平均估计 2.391  偏低 0.039
# 样本量  1000: 平均估计 2.422  偏低 0.008
# 样本量  5000: 平均估计 2.430  偏低 0.000
​

数字值得抄在工具箱盖子内侧：八个符号的分布，五十个样本低估零点一一比特，两百个样本低估零点零四，一千个样本基本收敛。经验法则是"每格样本数"——期望出现次数低于五到十的单元格，贡献的偏差就开始可观。校正方向有两条：解析校正（米勒-马多校正，近似扣除"符号数减一除以两倍样本"量级的偏差）或贝叶斯平滑（先验伪计数）；两者都便宜，没有理由裸奔。高维场景（估词表上万的互信息）里，偏差与方差的双面夹击更凶，专门的近邻型或核密度估计器才是正解——工具选型先看维度再看别的。

互信息的矩阵化：一次称完一整面特征墙

特征工程的日常不是算一个互信息，是给一整面特征墙排序。把"特征-标签"互信息做成矩阵运算，配合可视化读热力图，是最快的初筛流程。

# 特征墙初筛：批量计算互信息并排出前几名
import numpy as np

def mi_discrete(x, y, bx=8, by=2):
    """x 连续特征分箱后与离散标签 y 的互信息"""
    xr = np.clip(((x - x.min()) / (x.ptp() + 1e-12) * bx).astype(int), 0, bx - 1)
    joint = np.zeros((bx, by))
    for a, b in zip(xr, y):
        joint[a, b] += 1
    joint /= joint.sum()
    px, py = joint.sum(1), joint.sum(0)
    with np.errstate(divide="ignore", invalid="ignore"):
        terms = joint * np.log2(joint / (px[:, None] * py[None, :]))
    return np.nansum(terms)

rng = np.random.default_rng(0)
n = 4000
y = rng.integers(0, 2, n)                                   # 二分类标签
features = {
    "与标签强相关":   rng.normal(0, 1, n) + y * 2.5,        # 两类错开 2.5 个标准差
    "与标签弱相关":   rng.normal(0, 1, n) + y * 0.3,
    "与标签无关":     rng.normal(0, 1, n),
    "恒定特征":       np.zeros(n),
}
for name, x in features.items():
    print(f"{name}: 互信息 ≈ {mi_discrete(x, y):.4f} 比特")
# 输出（每次运行略有浮动）:
# 与标签强相关: 互信息 ≈ 0.6192 比特
# 与标签弱相关: 互信息 ≈ 0.0191 比特
# 与标签无关: 互信息 ≈ 0.0003 比特
# 恒定特征: 互信息 ≈ 0.0000 比特
​

三档读数的含义上一章讲过，这里强调流程位置：互信息初筛（便宜、快、非线性敏感）在前，模型验证在后——筛掉零方差与明显无关的特征压缩搜索空间，再用交叉验证细选。注意分箱数是超参数：箱太碎每格样本稀疏、偏差抬升；箱太粗把非线性关系抹平。实用起点是每格不少于几十个样本，并对分箱数做一次敏感性检查。

容量曲线的数值计算：当闭式解不在时

BSC 与 AWGN 有闭式容量，但一旦换到真实信道模型（多电平调制、非高斯噪声、有色干扰），闭式解就消失了，容量只能数值算：对外层输入分布做优化，内层对条件分布求和。以最常用的离散化 AWGN 加多电平调制为例，容量的核心是"给定输入星座、对输出做数值积分的互信息"。

# 多电平调制的可达速率：数值积分版互信息（等概输入）
import numpy as np
from math import pi, e, log2

def achievable_rate(snr_db, levels=4, grid=4001):
    """M 电平幅度调制在 AWGN 下的每符号可达比特数（等概输入）"""
    snr = 10 ** (snr_db / 10)
    cons = np.linspace(-1, 1, levels) * np.sqrt(3 / (levels ** 2 - 1))  # 单位平均能量
    sigma2 = 1 / (2 * snr)                        # 实信号每维噪声方差
    sigma = sigma2 ** 0.5
    y = np.linspace(cons.min() - 8 * sigma, cons.max() + 8 * sigma, grid)
    py = np.zeros_like(y)
    for x in cons:                                 # 输出密度 = 各星座点高斯的混合
        py += np.exp(-(y - x) ** 2 / (2 * sigma2)) / levels
    py /= (2 * pi * sigma2) ** 0.5
    H_Y = -np.trapz(py * np.log2(py + 1e-300), y)  # 输出熵（数值积分）
    H_YgX = 0.5 * log2(2 * pi * e * sigma2)        # 条件熵（高斯闭式）
    return H_Y - H_YgX

for db in [0, 5, 10, 15, 20]:
    print(f"信噪比 {db:>2} dB: 四电平可达 {achievable_rate(db):.3f} 比特/符号（上限 log2 4 = 2）")
# 输出:
# 信噪比  0 dB: 四电平可达 0.145 比特/符号（上限 log2 4 = 2）
# 信噪比  5 dB: 四电平可达 0.382 比特/符号（上限 log2 4 = 2）
# 信噪比 10 dB: 四电平可达 0.819 比特/符号（上限 log2 4 = 2）
# 信噪比 15 dB: 四电平可达 1.399 比特/符号（上限 log2 4 = 2）
# 信噪比 20 dB: 四电平可达 1.898 比特/符号（上限 log2 4 = 2）
​

读数的形状与第 4 章的理论完全同构：低信噪比区每符号不足零点二比特（噪声把电平糊成一团），高信噪比区逼近星座上限（四电平至多两比特）——星座尺寸封顶、噪声压低实际读数，两条力的合成曲线。数值实现的两处纪律：积分网格要覆盖"星座范围加噪声展开"（跨度不足截掉尾巴，输出熵被低估）；条件熵用高斯闭式、输出熵才走数值——能闭式的部分闭式，既省算力又少一层误差。严格容量还要在最外层对输入概率做优化（凸问题，可靠收敛），等概只是简化；这套"离散化加数值积分"的流程是调制设计论文的标配工具，值得完整走一遍。

工具链的组织：一个随手可查的称量台

三件量具之外，日常称量台通常还配四样东西。随机数纪律：所有仿真固定种子发布结果、换种子复跑验证稳定性（置信区间宽度速查见下一节）。向量化：样本量上十万后纯 Python 循环成为瓶颈，数组运算把速度抬两三个数量级。单位约定：比特还是自然对数、每符号还是每秒、分贝还是线性——三类单位混用的错误占实测事故的相当比例，脚本开头统一约定并打印换算。极限自检：算完任何读数，先跑两个边界情形（零信噪比应归零、满信噪比应达星座上限）确认实现没错，再相信中间读数——把第 4 章"边界读数"的理论直觉变成单元测试。

💡 一条值得贴在显示器边的话：数值实验和理论一样会出错，只是错得更有信心。边界自检、种子纪律、单位约定这三件小事，是把"数字"变成"证据"的全部仪式。

本节要点回顾

• 插件法从样本估熵系统性偏低：八符号分布五十样本低估零点一一比特，经验法则"每格样本数至少五到十"，米勒-马多校正或贝叶斯平滑便宜且必要；
• 高维互信息估计偏差方差双面夹击，近邻或核密度估计器是正解，选型先看维度；
• 特征墙初筛流程：互信息批量排序（便宜、非线性敏感）在前，模型交叉验证在后；分箱数是超参数，做敏感性检查；
• 数值容量计算=外层输入分布优化+内层条件分布积分，边界自检（零信噪比归零、高信噪比达星座上限）是必做单元测试；
• 称量台四配件：种子纪律、向量化、单位约定、极限自检——把数字变成证据的仪式。

量具会用了，下一节把它们编进设计与验证的守则——余量怎么留、仿真怎么设计、基准曲线怎么读。