尧图网页设计 
简介面向信息论与复杂系统分析需求这份传递熵计算资源为从事信号处理、通信网络及因果分析的研究人员与工程师提供了可直接使用的MATLAB实现内容紧密围绕传递熵、延迟时间与最大熵原理三大核心概念通过联合概率与条件熵计算帮助量化两个系统间信息的传递方向与程度。压缩包共3个文件以m源码为主并附带一张示意图说明算法流程整体仅20KB轻量便于快速部署代码覆盖数据读取、联合概率与条件熵计算、传递熵求解以及不同延迟时间下的扫描寻优能定位使传递熵最大化的时间间隔辅助研究者分析系统响应速度与耦合机制便于后续在此基础上做进一步算法改造与功能扩展。已有221人学习该资源读者可在此基础上结合自身实验数据开展验证进一步扩展数据预处理与可视化流程适用于通信链路效能评估、神经信号耦合分析、金融市场联动探测等多种实际任务。1. 传递熵计算是什么用信息流找“谁领先谁几拍”时序分析做到第三、四周很多人都会从“看相关性”转向“找因果性”。两个序列算出的皮尔逊相关系数再高也没法告诉你到底是 A 拽着 B还是 B 拖着 A滞后互相关能找到一个延迟峰值但面对非线性耦合和噪声干扰时这个峰经常出现在错误的方向上。传递熵计算解决的就是这个问题用条件互信息衡量“在已知接收方过去的前提下发送方的过去对接收方未来额外消除了多少不确定”从而给出一个有方向的、可扫描延迟的量。我第一次接触这套方案是在做新能源功率序列分析时天气系统、风功率和光伏出力互相纠缠线性方法全都在贡献伪相关最后是靠传递熵把主导方向刻出来的。这类项目经常以“传递熵计算”脚本包的形式出现压缩包里通常就是三块东西概率密度估计函数、延迟时间扫描逻辑、以及最大熵原理相关的辅助估计模块。阅读这篇文章的人不需要懂复杂数学只要会 Python 和 numpy跟着代码走完一遍就能在自己的双变量时序上算出 TE 曲线再配合最大熵谱估计缩小滞后搜索范围。下面我会从公式、实现、参数到坑位把整个流程拆给你看。2. 传递熵的公式拆解条件互信息与两个关键参数2.1 条件互信息视角下的传递熵传递熵 $TE_{Y \to X}$ 的定义用条件互信息写最直观$$ TE_{Y \to X}(s) I(X_{ns}; Y_n \mid X_n) \sum p(x_{ns}, x_n, y_n) \log \frac{p(x_{ns} \mid x_n, y_n)}{p(x_{ns} \mid x_n)} $$这里 $Y$ 是发送方$X$ 是接收方$s$ 是延迟步数。右侧的条件概率比里分子是同时知道 $X_n$ 和 $Y_n$ 后对 $X_{ns}$ 的预测能力分母是只知道 $X_n$ 时的预测能力。如果 $Y$ 对 $X$ 的将来没有任何额外信息条件概率相等TE 为零一旦存在耦合信息差就会让 TE 变成正数。相比互信息传递熵天然是非对称的。$TE_{Y \to X}$ 和 $TE_{X \to Y}$ 可能在数值上有明显差异这正是判断因果方向的依据。需要注意这种“因果”是预测意义上的不是实验干预意义上的因果所以更严谨的说法是“发现方向性的信息流而不是证明物理因果”。2.2 延迟时间 s 决定信息传递的时滞延迟时间 $s$ 是传递熵计算里最重要的参数。$s$ 的意义是发送方今天状态贡献的信息需要经过多少个采样间隔才在接收方未来状态中体现出来。实际使用中我没有只算一个 s而是扫描一个范围例如s 1..20画出 $TE_{Y \to X}(s)$ 的曲线。曲线出现明显峰值的那个位置就是估计的传递延迟。如果你有嵌入维数概念这里还要注意一点公式里的 $X_n$ 和 $Y_n$ 往往只是一个标量可真实物理系统可能是高阶动力学。一种扩展做法是把 $X_n$ 替换成过去的嵌入向量 $X_n^{(m)}$把 $Y_n$ 替换成 $Y_n^{(l)}$得到更全面的条件互信息。代价是联合概率估计的维度上升数据量需求急剧增大所以工程上经常先跑一阶近似把峰值方向选出来再针对候选延迟做高阶验证。2.3 平稳化与去趋势预处理必须在前传递熵本质上是一堆概率密度的比值它对联合分布的平稳性很敏感。趋势项会让 $x_n$ 和 $x_{ns}$ 的均值位置漂移联合密度结构被时间趋势主导算出来的 TE 往往要么虚高要么完全掩盖真实耦合。我一般按三步走第一步对每个序列做一阶差分消除线性趋势第二步如果还有明显的周期性用滚动窗口减去局部均值第三步把两个序列都减均值、除以标准差让数值尺度统一。需要注意差分会放大高频噪声如果信号本身信噪比很低可以考虑换成符号化方式用排序模式代替原始数值这个我在最后一章再展开。3. 用原生 Python 把传递熵计算跑通3.1 直方图最小实现三维联合密度估计最简单的传递熵实现是直方图法。我们把原始序列拆成三列向量接收方过去 $x_n$、发送方过去 $y_n$、接收方未来 $x_{ns}$然后估计三维联合概率分布再通过降维得到条件概率比。下面是完整可跑的最小实现。import numpy as np def hist_te(x, y, lag, bins8): x np.asarray(x, dtypefloat).ravel() y np.asarray(y, dtypefloat).ravel() n len(x) - lag if n 0: raise ValueError(flag{lag} 超过了序列长度) x_past x[:n] y_past y[:n] x_fut x[lag:] # 三个向量长度都是 n samples np.stack([x_past, y_past, x_fut], axis1) # 三维直方图维度 0x_past, 1y_past, 2x_fut counts, edges np.histogramdd(samples, binsbins) P3 counts / counts.sum() # 二维边际概率 Pxy P3.sum(axis2) # p(x_past, y_past) Pzx P3.sum(axis1) # p(x_past, x_fut) Px P3.sum(axis(1, 2)) # p(x_past) # p(x_fut | x_past, y_past) cond_y P3 / Pxy[:, :, None] # p(x_fut | x_past) cond_x Pzx / Px[:, None] # 条件概率比NaN/inf 留在 P3 0 的位置之后被忽略 with np.errstate(divideignore, invalidignore): ratio cond_y / cond_x[:, None, :] log_ratio np.log(ratio) te np.sum(P3 * log_ratio, where(P3 0)) return te这段代码里最需要注意的是np.histogramdd返回的轴顺序。传入samples的三列分别是x_past、y_past、x_fut所以三维数组第 0 维是x_past的格子第 1 维是y_past的格子第 2 维是x_fut的格子。降维求和时axis2消掉未来维度得到发送方和接收方联合axis1消掉发送方得到接收方过去与未来联合axis(1,2)消掉发送方和未来得到接收方过去的边际概率。计算ratio之前我特意用了np.errstate因为在稀疏直方图里cond_x可能为零产生除零警告。真正参与求和时where(P3 0)保证只有联合概率非零的格子才计入 TE避免了0 * inf产生 NaN。3.2 KDE 版本用核密度估计替代格子直方图的主要问题是格子边界不连续高维情况下样本稀疏边际概率容易塌缩。KDE 的解法是把事件点周围的信息用核函数扩散开密度估计更光滑对中等长度数据更稳。sklearn 的 KernelDensity 可以直接用来估计对数密度简洁地算出 TE。from sklearn.neighbors import KernelDensity def kde_te(x, y, lag1, bandwidth0.3): x np.asarray(x, dtypefloat).ravel() y np.asarray(y, dtypefloat).ravel() n len(x) - lag x_past x[:n] y_past y[:n] x_fut x[lag:] X3 np.c_[x_past, y_past, x_fut] X_xy np.c_[x_past, y_past] X_xz np.c_[x_past, x_fut] X_x x_past.reshape(-1, 1) kd3 KernelDensity(bandwidthbandwidth).fit(X3) kd_xy KernelDensity(bandwidthbandwidth).fit(X_xy) kd_xz KernelDensity(bandwidthbandwidth).fit(X_xz) kd_x KernelDensity(bandwidthbandwidth).fit(X_x) # 在样本点本身处评估 log 密度 log3 kd3.score_samples(X3) log_xy kd_xy.score_samples(X_xy) log_xz kd_xz.score_samples(X_xz) log_x kd_x.score_samples(X_x) # TE E[ log p(x_fut,x_past,y_past) log p(x_past) # - log p(x_past,y_past) - log p(x_past,x_fut) ] te np.mean(log3 log_x - log_xy - log_xz) return te这个版本巧妙绕开了显式计算条件概率直接在样本点上用对数密度做期望。为什么可以这样因为条件概率比展开后等于联合密度乘积比所以 TE 近似为四个 KDE 对数密度值的差值平均。KDE 在高斯核下对数据尺度敏感所以在调用前最好把序列标准化。坏处也不是没有KDE 在训练样本点本身做评估会存在轻微的正偏差因为每个点对自己有贡献。对比两个方向的 TE 时这个偏差大致对称影响不大如果要做严格的显著性检验最好把数据集分成训练和评估两半但那样样本量要翻倍短序列容易撑不住。3.3 参数表bins、bandwidth、样本量与计算量参数直方图版本KDE 版本我的默认经验离散化参数binsbandwidthbins 取 8~12bandwidth 取 0.2~0.5 或 Scott 规则最小样本量500 起步800 起步希望更稳选 2000少于 300 不要做三维直方图计算复杂度随 bins 的立方增长训练 O(N log N)评估 O(N)N5000 时直方图秒级KDE 约分钟级敏感性对边界位置敏感对带宽极其敏感参考下一章的等概率分箱表格里最需要记住的是直方图不只是看格子数还要看每个格子里的期望样本数。如果bins20但数据只有 1000 点三维空间就有 8000 个格子平均一个格子只有 0.125 个样本结果全是噪声。KDE 的带宽同理带宽太小密度函数变成一堆尖刺TE 完全由个体差异主导带宽太大所有概率分布都变得接近均匀传递熵信息被抹平。3.4 用合成信号做自检先确认方向没有被代码搞反任何新写的 TE 函数都要先跑合成数据验证方向。我常用的生成方式是构造一个严格的单向耦合让 $x$ 是白噪声$y$ 在延迟两步之后跟随 $x$但反向 $x$ 不跟随 $y$。这样理论上 $TE_{X \to Y}$ 应该明显大于 $TE_{Y \to X}$。rng np.random.default_rng(42) n 3000 x rng.standard_normal(n) y rng.standard_normal(n) for t in range(2, n): y[t] 0.6 * x[t - 2] 0.2 * y[t - 1] 0.1 * rng.standard_normal() te_x_to_y hist_te(x, y, lag2, bins8) te_y_to_x hist_te(y, x, lag2, bins8) print(fX - Y: {te_x_to_y:.4f}) print(fY - X: {te_y_to_x:.4f})如果代码写错轴顺序最常见的结果是两个方向几乎相等或者反向比正向还大。直方图法在小样本下总是会高估 TE所以不必追求接近理论零值只要两个方向有显著差异且方向正确就可以继续往下做延迟扫描。4. 延迟时间扫描与最大熵原理把短样本问题交给约束优化4.1 别只看单点 TE扫一段滞后曲线真实耦合不是只在一个整数延迟上存在信号经过调制、传播和滤波后信息可能分布在连续几个滞后上。我一般是取lag_max 20或者 1/4 的采样数量从 1 扫到lag_max得到一个 TE 曲线。曲线峰值所在的位置最有价值但峰的形状同样重要如果曲线只在某个延迟处尖锐凸起说明动态关系明确如果整个曲线都高说明可能是共驱动因素需要做条件传递熵。为什么要先扫延迟而不是直接猜一个因为传递熵对滞后很敏感。把 s 设小可能没覆盖到真实传播时延把 s 设大条件概率中接收方过去的信息已经包含了部分发送方影响TE 会被低估。只有在正确的 s 上接收方过去不会“抢走”发送方未来信息TE 的额外信息量才是最大的。4.2 最大熵原理是什么为什么适合短序列谱估计延迟扫描范围不是越宽越好。如果信号里有明显周期成分最自然的做法是把最大扫描延迟设成主周期长度的一倍到两倍避免在无关滞后上反复试。问题是短序列的傅里叶周期图旁瓣大谱峰不稳定这时候就轮到最大熵原理上场。最大熵原理的基本思想是在已知约束条件下选择熵最大的概率分布作为估计结果因为这是对未知部分最“审慎”的假设。应用到时间序列里假设已知自相关函数前 p 个值其他一切未知那么熵最大化得到的谱正是 AR(p) 模型对应的功率谱。这也是为什么最大熵谱估计在短样本上往往比经典周期图更平滑、峰更稳定——它不假设观测窗口之外数据为零而是用满足约束的最随机延拓来推断。4.3 用 Yule-Walker 方程实现最大熵谱估计并圈定延迟范围一段能直接跑的 Yule-Walker 最大熵谱估计函数如下from scipy.linalg import toeplitz def maxent_spectrum(x, order20, fs1.0, n_freqs256): x np.asarray(x, dtypefloat).ravel() x x - np.mean(x) N len(x) # 自相关函数除以 N 得到无偏估计的近似 r np.correlate(x, x, modefull)[N - 1:] / N # Yule-Walker 方程R * a -r[1:] R toeplitz(r[:order]) rhs -r[1:order 1] ar_coef np.linalg.solve(R, rhs) # AR 系数第一个是 1 a np.r_[1.0, ar_coef] freqs np.linspace(0, fs / 2, n_freqs) A np.zeros_like(freqs, dtypecomplex) for k, coef in enumerate(a): A coef * np.exp(-2j * np.pi * k * freqs / fs) # 最大熵谱密度与 |A(freq)|^2 成反比 spec 1.0 / (np.abs(A) ** 2) return freqs, spec # 使用示例 freqs, spec maxent_spectrum(x, order30, fs1.0) peak_freq freqs[np.argmax(spec)] print(f主导频率: {peak_freq:.4f} Hz, 对应周期约 {1/peak_freq:.1f} 个采样点)参数order是 AR 阶数控制谱的复杂度。order 太小谱峰太平看不出周期order 太大会把噪声当成真实谱峰经验值是order N / 10我习惯从 20 起步。得到主周期 T 之后把延迟扫描范围设置成 1 到 T 之间最多到 2*T就能避开大量无效滞后。这套流程特别适合短样本数据最大熵谱估计在 500 点时依然能给出可用的主峰而周期图可能已经出现多个不相关旁瓣。注意它给出的是发送方或接收方各自内部的主导周期不是交互延迟交互延迟还是要靠 TE 曲线峰定位谱估计只负责缩小候选区。4.4 用最大熵思想设置三维直方图的等概率边界很多人在 hist_te 里传bins8默认等宽划分但等宽对数据分布不均匀的场景很吃亏。一段信号如果大部分集中在 0 附近等宽格子里外围格子完全为空三维联合概率的估计效率变差。最大熵思想在这里给出一个实用改良让每个维度上的边界点是样本分位数使每个一维区间包含大致相等的样本数。def quantile_edges(data, n_bins): q np.linspace(0, 1, n_bins 1) return np.quantile(data, q) # 替换 hist_te 中的 hist_te 调用方法 n_bins 8 ex quantile_edges(x_past, n_bins) ey quantile_edges(y_past, n_bins) ez quantile_edges(x_fut, n_bins) counts, _ np.histogramdd(samples, bins[ex, ey, ez])这里没有真正去求解带约束的熵最大化但等概率分箱在离散化过程中最小化信息损失与最大熵原则的方向一致。相比等宽分箱它对尾部样本更友好能够保住少量极端事件对信息流的贡献。注意x_past和x_fut虽然来自同一个序列但对应时间段不同经验上建议各自算边界不要共用一份分位数直接用同一份边界也不会错只是边缘区间可能样本分配不均。5. 传递熵计算避坑指南五个常见翻车场景5.1 零滞后伪峰算出来的 TE 在 s0 处特别大现象把lag0放进 TE 函数得到的数值比所有正延迟都大而且两个方向几乎对称。原因零滞后时$x_{ns}$ 就是 $x_n$条件互信息实际上在测“同一时刻两个变量的同步耦合”而不是预测意义下的信息传递。如果两个序列受同一个外部驱动比如同一段电网电压波动或者同一个天气过程零滞后同步会被误判成双向信息流。解决永远从lag1开始扫描。如果确实需要评估瞬时同步把结果单独标注为“同步强度”不要和传递熵混为一谈。另一个补救是用置换检验把发送方序列随机循环平移后重算 TE过滤掉共同趋势带来的假阳性。5.2 样本长度不够直方图结果像随机数现象用 200 个点跑完以后TE 曲线在多个延迟上来回跳调整 bins 数值结论彻底反转。原因三维直方图要估计的联合概率空间很大。200 个样本放进 8^3 512 个格子里平均每个格子不到 1 个样本统计波动极大bins 一变格子边界一变密度结构立刻就变。解决我的底线是 N 500 起步再多也不嫌多。如果只有 200 点优先用 KDE 或者符号化 TE不要硬上三维直方图。另一个折中是只估计二维条件互信息减少一个概率维度但那样会牺牲高阶耦合信息。5.3 bin 数固定导致稀疏矩阵和除零现象代码报RuntimeWarning: divide by zero最后 TE 结果是 NaN 或者无穷大。原因固定 bins 下很多格子的counts是零cond_x分母出现零。虽然用了np.errstate屏蔽警告但如果在处理前给P3加了1e-12平滑会把所有空格子都当成极小概率参与 log 计算噪声被线性放大。解决不要在密度估计上人为加小常数而是用where(P3 0)掩码处理对数项。如果要用平滑也要用 Dirichlet 先验密度比如在counts alpha之后重新归一化alpha 取 0.5 这类值而不是随意加一个常数。5.4 带宽选择让 KDE 方向性消失现象KDE 版本算出来两个方向的 TE 几乎一样或者数值大到离谱。原因bandwidth 设得太小时KDE 密度函数退化成一系列脉冲每个样本点只看得到自己条件概率比趋近于 1TE 被高估bandwidth 设得太大时所有密度函数都变成同一个扁平高斯条件概率比趋近于 1TE 又被低估。两个极端都会抹平方向差异。解决用 Scott 规则给一个起点bandwidth N ** (-1 / (d 4))其中 d 是维度。三维联合密度用 d3二维边际用 d2但 KDE 函数里我用的是同一个带宽简化处理。靠谱的做法是分别用 0.5、1、2 倍默认带宽各跑一遍看 TE 峰值位置是否稳定峰值位置能稳定数值有差别结论才值得信任。5.5 不做置换检验任何正数 TE 都没有意义现象算出的 TE 0.02看起来是正数于是宣布找到了方向性耦合。原因传递熵的估计量是正偏的直方图法和 KDE 法在有限样本下都会产生正的“基础噪声”。两个白噪声序列也能算出一个明显大于零的 TE这并不代表有信息流。解决做零假设检验。把发送方序列循环平移随机步破坏原始配对关系再重新计算 TE重复 200 到 1000 次得到零分布。如果真实 TE 超过了零分布的 95% 分位点才可以说存在显著信息流。循环平移比简单打乱好因为它能保留发送方序列内部的时间相关性避免因为自相关结构改变而产生假结论。6. 进阶用法多变量条件 TE 与稳定性验证6.1 条件传递熵剔除公共驱动当第三个序列 Z 同时驱动 X 和 Y 时普通 TE 会把虚假方向算出来。条件传递熵在估计条件概率时把 Z 的过去也加入条件集合也就是把p(x_{ns} | x_n, y_n)改成p(x_{ns} | x_n, y_n, z_n)并与不包含 y_n 的条件概率做比值。实现上就是在直方图或 KDE 里多加一列数据样本需求会大幅上升所以高阶场景建议直接用符号化。6.2 符号化排序近似符号化 TE 不估计连续密度而是把三维样本按数值排序。对每个时刻 t看x_past、y_past、x_fut的大小顺序用排列模式代替原始值然后在离散模式上计数计算条件互信息。这个方法对异常值和噪声极稳参数只有一个延迟 s不需要调带宽。缺点是三变量的排序模式最多 6 种信息分辨率有限适合作为 TE 结果的交叉验证。6.3 稳定性验证我自己的三条结局检查我现在拿到任何 TE 结果都要求自己把这三件事做齐第一画 TE 随延迟变化的曲线而不是只报一个峰值第二把数据切成两段分别计算看峰值延迟是否一致第三至少做 500 次置换检验把原始 TE 和零分布画在同一张图上。如果峰值移动超过两个采样点或者 p 值大于 0.05我宁可不出报告。多年下来这个习惯帮我抓住了很多本来会写进结论里的假因果。实现层面还有一个实用技巧把延迟扫描和带宽测试封装成一个循环输出一个三列数组包含延迟、两个方向的 TE、以及置换 p 值这样每个项目只要跑一次结果就能直接进周报和幻灯片。做传递熵计算时“算得出来”永远不是终点“稳定可复现”才是。希望帮到你。本文还有配套的精品资源点击获取 传递熵计算：用条件互信息识别时序因果方向
简介面向信息论与复杂系统分析需求这份传递熵计算资源为从事信号处理、通信网络及因果分析的研究人员与工程师提供了可直接使用的MATLAB实现内容紧密围绕传递熵、延迟时间与最大熵原理三大核心概念通过联合概率与条件熵计算帮助量化两个系统间信息的传递方向与程度。压缩包共3个文件以m源码为主并附带一张示意图说明算法流程整体仅20KB轻量便于快速部署代码覆盖数据读取、联合概率与条件熵计算、传递熵求解以及不同延迟时间下的扫描寻优能定位使传递熵最大化的时间间隔辅助研究者分析系统响应速度与耦合机制便于后续在此基础上做进一步算法改造与功能扩展。已有221人学习该资源读者可在此基础上结合自身实验数据开展验证进一步扩展数据预处理与可视化流程适用于通信链路效能评估、神经信号耦合分析、金融市场联动探测等多种实际任务。1. 传递熵计算是什么用信息流找“谁领先谁几拍”时序分析做到第三、四周很多人都会从“看相关性”转向“找因果性”。两个序列算出的皮尔逊相关系数再高也没法告诉你到底是 A 拽着 B还是 B 拖着 A滞后互相关能找到一个延迟峰值但面对非线性耦合和噪声干扰时这个峰经常出现在错误的方向上。传递熵计算解决的就是这个问题用条件互信息衡量“在已知接收方过去的前提下发送方的过去对接收方未来额外消除了多少不确定”从而给出一个有方向的、可扫描延迟的量。我第一次接触这套方案是在做新能源功率序列分析时天气系统、风功率和光伏出力互相纠缠线性方法全都在贡献伪相关最后是靠传递熵把主导方向刻出来的。这类项目经常以“传递熵计算”脚本包的形式出现压缩包里通常就是三块东西概率密度估计函数、延迟时间扫描逻辑、以及最大熵原理相关的辅助估计模块。阅读这篇文章的人不需要懂复杂数学只要会 Python 和 numpy跟着代码走完一遍就能在自己的双变量时序上算出 TE 曲线再配合最大熵谱估计缩小滞后搜索范围。下面我会从公式、实现、参数到坑位把整个流程拆给你看。2. 传递熵的公式拆解条件互信息与两个关键参数2.1 条件互信息视角下的传递熵传递熵 $TE_{Y \to X}$ 的定义用条件互信息写最直观$$ TE_{Y \to X}(s) I(X_{ns}; Y_n \mid X_n) \sum p(x_{ns}, x_n, y_n) \log \frac{p(x_{ns} \mid x_n, y_n)}{p(x_{ns} \mid x_n)} $$这里 $Y$ 是发送方$X$ 是接收方$s$ 是延迟步数。右侧的条件概率比里分子是同时知道 $X_n$ 和 $Y_n$ 后对 $X_{ns}$ 的预测能力分母是只知道 $X_n$ 时的预测能力。如果 $Y$ 对 $X$ 的将来没有任何额外信息条件概率相等TE 为零一旦存在耦合信息差就会让 TE 变成正数。相比互信息传递熵天然是非对称的。$TE_{Y \to X}$ 和 $TE_{X \to Y}$ 可能在数值上有明显差异这正是判断因果方向的依据。需要注意这种“因果”是预测意义上的不是实验干预意义上的因果所以更严谨的说法是“发现方向性的信息流而不是证明物理因果”。2.2 延迟时间 s 决定信息传递的时滞延迟时间 $s$ 是传递熵计算里最重要的参数。$s$ 的意义是发送方今天状态贡献的信息需要经过多少个采样间隔才在接收方未来状态中体现出来。实际使用中我没有只算一个 s而是扫描一个范围例如s 1..20画出 $TE_{Y \to X}(s)$ 的曲线。曲线出现明显峰值的那个位置就是估计的传递延迟。如果你有嵌入维数概念这里还要注意一点公式里的 $X_n$ 和 $Y_n$ 往往只是一个标量可真实物理系统可能是高阶动力学。一种扩展做法是把 $X_n$ 替换成过去的嵌入向量 $X_n^{(m)}$把 $Y_n$ 替换成 $Y_n^{(l)}$得到更全面的条件互信息。代价是联合概率估计的维度上升数据量需求急剧增大所以工程上经常先跑一阶近似把峰值方向选出来再针对候选延迟做高阶验证。2.3 平稳化与去趋势预处理必须在前传递熵本质上是一堆概率密度的比值它对联合分布的平稳性很敏感。趋势项会让 $x_n$ 和 $x_{ns}$ 的均值位置漂移联合密度结构被时间趋势主导算出来的 TE 往往要么虚高要么完全掩盖真实耦合。我一般按三步走第一步对每个序列做一阶差分消除线性趋势第二步如果还有明显的周期性用滚动窗口减去局部均值第三步把两个序列都减均值、除以标准差让数值尺度统一。需要注意差分会放大高频噪声如果信号本身信噪比很低可以考虑换成符号化方式用排序模式代替原始数值这个我在最后一章再展开。3. 用原生 Python 把传递熵计算跑通3.1 直方图最小实现三维联合密度估计最简单的传递熵实现是直方图法。我们把原始序列拆成三列向量接收方过去 $x_n$、发送方过去 $y_n$、接收方未来 $x_{ns}$然后估计三维联合概率分布再通过降维得到条件概率比。下面是完整可跑的最小实现。import numpy as np def hist_te(x, y, lag, bins8): x np.asarray(x, dtypefloat).ravel() y np.asarray(y, dtypefloat).ravel() n len(x) - lag if n 0: raise ValueError(flag{lag} 超过了序列长度) x_past x[:n] y_past y[:n] x_fut x[lag:] # 三个向量长度都是 n samples np.stack([x_past, y_past, x_fut], axis1) # 三维直方图维度 0x_past, 1y_past, 2x_fut counts, edges np.histogramdd(samples, binsbins) P3 counts / counts.sum() # 二维边际概率 Pxy P3.sum(axis2) # p(x_past, y_past) Pzx P3.sum(axis1) # p(x_past, x_fut) Px P3.sum(axis(1, 2)) # p(x_past) # p(x_fut | x_past, y_past) cond_y P3 / Pxy[:, :, None] # p(x_fut | x_past) cond_x Pzx / Px[:, None] # 条件概率比NaN/inf 留在 P3 0 的位置之后被忽略 with np.errstate(divideignore, invalidignore): ratio cond_y / cond_x[:, None, :] log_ratio np.log(ratio) te np.sum(P3 * log_ratio, where(P3 0)) return te这段代码里最需要注意的是np.histogramdd返回的轴顺序。传入samples的三列分别是x_past、y_past、x_fut所以三维数组第 0 维是x_past的格子第 1 维是y_past的格子第 2 维是x_fut的格子。降维求和时axis2消掉未来维度得到发送方和接收方联合axis1消掉发送方得到接收方过去与未来联合axis(1,2)消掉发送方和未来得到接收方过去的边际概率。计算ratio之前我特意用了np.errstate因为在稀疏直方图里cond_x可能为零产生除零警告。真正参与求和时where(P3 0)保证只有联合概率非零的格子才计入 TE避免了0 * inf产生 NaN。3.2 KDE 版本用核密度估计替代格子直方图的主要问题是格子边界不连续高维情况下样本稀疏边际概率容易塌缩。KDE 的解法是把事件点周围的信息用核函数扩散开密度估计更光滑对中等长度数据更稳。sklearn 的 KernelDensity 可以直接用来估计对数密度简洁地算出 TE。from sklearn.neighbors import KernelDensity def kde_te(x, y, lag1, bandwidth0.3): x np.asarray(x, dtypefloat).ravel() y np.asarray(y, dtypefloat).ravel() n len(x) - lag x_past x[:n] y_past y[:n] x_fut x[lag:] X3 np.c_[x_past, y_past, x_fut] X_xy np.c_[x_past, y_past] X_xz np.c_[x_past, x_fut] X_x x_past.reshape(-1, 1) kd3 KernelDensity(bandwidthbandwidth).fit(X3) kd_xy KernelDensity(bandwidthbandwidth).fit(X_xy) kd_xz KernelDensity(bandwidthbandwidth).fit(X_xz) kd_x KernelDensity(bandwidthbandwidth).fit(X_x) # 在样本点本身处评估 log 密度 log3 kd3.score_samples(X3) log_xy kd_xy.score_samples(X_xy) log_xz kd_xz.score_samples(X_xz) log_x kd_x.score_samples(X_x) # TE E[ log p(x_fut,x_past,y_past) log p(x_past) # - log p(x_past,y_past) - log p(x_past,x_fut) ] te np.mean(log3 log_x - log_xy - log_xz) return te这个版本巧妙绕开了显式计算条件概率直接在样本点上用对数密度做期望。为什么可以这样因为条件概率比展开后等于联合密度乘积比所以 TE 近似为四个 KDE 对数密度值的差值平均。KDE 在高斯核下对数据尺度敏感所以在调用前最好把序列标准化。坏处也不是没有KDE 在训练样本点本身做评估会存在轻微的正偏差因为每个点对自己有贡献。对比两个方向的 TE 时这个偏差大致对称影响不大如果要做严格的显著性检验最好把数据集分成训练和评估两半但那样样本量要翻倍短序列容易撑不住。3.3 参数表bins、bandwidth、样本量与计算量参数直方图版本KDE 版本我的默认经验离散化参数binsbandwidthbins 取 8~12bandwidth 取 0.2~0.5 或 Scott 规则最小样本量500 起步800 起步希望更稳选 2000少于 300 不要做三维直方图计算复杂度随 bins 的立方增长训练 O(N log N)评估 O(N)N5000 时直方图秒级KDE 约分钟级敏感性对边界位置敏感对带宽极其敏感参考下一章的等概率分箱表格里最需要记住的是直方图不只是看格子数还要看每个格子里的期望样本数。如果bins20但数据只有 1000 点三维空间就有 8000 个格子平均一个格子只有 0.125 个样本结果全是噪声。KDE 的带宽同理带宽太小密度函数变成一堆尖刺TE 完全由个体差异主导带宽太大所有概率分布都变得接近均匀传递熵信息被抹平。3.4 用合成信号做自检先确认方向没有被代码搞反任何新写的 TE 函数都要先跑合成数据验证方向。我常用的生成方式是构造一个严格的单向耦合让 $x$ 是白噪声$y$ 在延迟两步之后跟随 $x$但反向 $x$ 不跟随 $y$。这样理论上 $TE_{X \to Y}$ 应该明显大于 $TE_{Y \to X}$。rng np.random.default_rng(42) n 3000 x rng.standard_normal(n) y rng.standard_normal(n) for t in range(2, n): y[t] 0.6 * x[t - 2] 0.2 * y[t - 1] 0.1 * rng.standard_normal() te_x_to_y hist_te(x, y, lag2, bins8) te_y_to_x hist_te(y, x, lag2, bins8) print(fX - Y: {te_x_to_y:.4f}) print(fY - X: {te_y_to_x:.4f})如果代码写错轴顺序最常见的结果是两个方向几乎相等或者反向比正向还大。直方图法在小样本下总是会高估 TE所以不必追求接近理论零值只要两个方向有显著差异且方向正确就可以继续往下做延迟扫描。4. 延迟时间扫描与最大熵原理把短样本问题交给约束优化4.1 别只看单点 TE扫一段滞后曲线真实耦合不是只在一个整数延迟上存在信号经过调制、传播和滤波后信息可能分布在连续几个滞后上。我一般是取lag_max 20或者 1/4 的采样数量从 1 扫到lag_max得到一个 TE 曲线。曲线峰值所在的位置最有价值但峰的形状同样重要如果曲线只在某个延迟处尖锐凸起说明动态关系明确如果整个曲线都高说明可能是共驱动因素需要做条件传递熵。为什么要先扫延迟而不是直接猜一个因为传递熵对滞后很敏感。把 s 设小可能没覆盖到真实传播时延把 s 设大条件概率中接收方过去的信息已经包含了部分发送方影响TE 会被低估。只有在正确的 s 上接收方过去不会“抢走”发送方未来信息TE 的额外信息量才是最大的。4.2 最大熵原理是什么为什么适合短序列谱估计延迟扫描范围不是越宽越好。如果信号里有明显周期成分最自然的做法是把最大扫描延迟设成主周期长度的一倍到两倍避免在无关滞后上反复试。问题是短序列的傅里叶周期图旁瓣大谱峰不稳定这时候就轮到最大熵原理上场。最大熵原理的基本思想是在已知约束条件下选择熵最大的概率分布作为估计结果因为这是对未知部分最“审慎”的假设。应用到时间序列里假设已知自相关函数前 p 个值其他一切未知那么熵最大化得到的谱正是 AR(p) 模型对应的功率谱。这也是为什么最大熵谱估计在短样本上往往比经典周期图更平滑、峰更稳定——它不假设观测窗口之外数据为零而是用满足约束的最随机延拓来推断。4.3 用 Yule-Walker 方程实现最大熵谱估计并圈定延迟范围一段能直接跑的 Yule-Walker 最大熵谱估计函数如下from scipy.linalg import toeplitz def maxent_spectrum(x, order20, fs1.0, n_freqs256): x np.asarray(x, dtypefloat).ravel() x x - np.mean(x) N len(x) # 自相关函数除以 N 得到无偏估计的近似 r np.correlate(x, x, modefull)[N - 1:] / N # Yule-Walker 方程R * a -r[1:] R toeplitz(r[:order]) rhs -r[1:order 1] ar_coef np.linalg.solve(R, rhs) # AR 系数第一个是 1 a np.r_[1.0, ar_coef] freqs np.linspace(0, fs / 2, n_freqs) A np.zeros_like(freqs, dtypecomplex) for k, coef in enumerate(a): A coef * np.exp(-2j * np.pi * k * freqs / fs) # 最大熵谱密度与 |A(freq)|^2 成反比 spec 1.0 / (np.abs(A) ** 2) return freqs, spec # 使用示例 freqs, spec maxent_spectrum(x, order30, fs1.0) peak_freq freqs[np.argmax(spec)] print(f主导频率: {peak_freq:.4f} Hz, 对应周期约 {1/peak_freq:.1f} 个采样点)参数order是 AR 阶数控制谱的复杂度。order 太小谱峰太平看不出周期order 太大会把噪声当成真实谱峰经验值是order N / 10我习惯从 20 起步。得到主周期 T 之后把延迟扫描范围设置成 1 到 T 之间最多到 2*T就能避开大量无效滞后。这套流程特别适合短样本数据最大熵谱估计在 500 点时依然能给出可用的主峰而周期图可能已经出现多个不相关旁瓣。注意它给出的是发送方或接收方各自内部的主导周期不是交互延迟交互延迟还是要靠 TE 曲线峰定位谱估计只负责缩小候选区。4.4 用最大熵思想设置三维直方图的等概率边界很多人在 hist_te 里传bins8默认等宽划分但等宽对数据分布不均匀的场景很吃亏。一段信号如果大部分集中在 0 附近等宽格子里外围格子完全为空三维联合概率的估计效率变差。最大熵思想在这里给出一个实用改良让每个维度上的边界点是样本分位数使每个一维区间包含大致相等的样本数。def quantile_edges(data, n_bins): q np.linspace(0, 1, n_bins 1) return np.quantile(data, q) # 替换 hist_te 中的 hist_te 调用方法 n_bins 8 ex quantile_edges(x_past, n_bins) ey quantile_edges(y_past, n_bins) ez quantile_edges(x_fut, n_bins) counts, _ np.histogramdd(samples, bins[ex, ey, ez])这里没有真正去求解带约束的熵最大化但等概率分箱在离散化过程中最小化信息损失与最大熵原则的方向一致。相比等宽分箱它对尾部样本更友好能够保住少量极端事件对信息流的贡献。注意x_past和x_fut虽然来自同一个序列但对应时间段不同经验上建议各自算边界不要共用一份分位数直接用同一份边界也不会错只是边缘区间可能样本分配不均。5. 传递熵计算避坑指南五个常见翻车场景5.1 零滞后伪峰算出来的 TE 在 s0 处特别大现象把lag0放进 TE 函数得到的数值比所有正延迟都大而且两个方向几乎对称。原因零滞后时$x_{ns}$ 就是 $x_n$条件互信息实际上在测“同一时刻两个变量的同步耦合”而不是预测意义下的信息传递。如果两个序列受同一个外部驱动比如同一段电网电压波动或者同一个天气过程零滞后同步会被误判成双向信息流。解决永远从lag1开始扫描。如果确实需要评估瞬时同步把结果单独标注为“同步强度”不要和传递熵混为一谈。另一个补救是用置换检验把发送方序列随机循环平移后重算 TE过滤掉共同趋势带来的假阳性。5.2 样本长度不够直方图结果像随机数现象用 200 个点跑完以后TE 曲线在多个延迟上来回跳调整 bins 数值结论彻底反转。原因三维直方图要估计的联合概率空间很大。200 个样本放进 8^3 512 个格子里平均每个格子不到 1 个样本统计波动极大bins 一变格子边界一变密度结构立刻就变。解决我的底线是 N 500 起步再多也不嫌多。如果只有 200 点优先用 KDE 或者符号化 TE不要硬上三维直方图。另一个折中是只估计二维条件互信息减少一个概率维度但那样会牺牲高阶耦合信息。5.3 bin 数固定导致稀疏矩阵和除零现象代码报RuntimeWarning: divide by zero最后 TE 结果是 NaN 或者无穷大。原因固定 bins 下很多格子的counts是零cond_x分母出现零。虽然用了np.errstate屏蔽警告但如果在处理前给P3加了1e-12平滑会把所有空格子都当成极小概率参与 log 计算噪声被线性放大。解决不要在密度估计上人为加小常数而是用where(P3 0)掩码处理对数项。如果要用平滑也要用 Dirichlet 先验密度比如在counts alpha之后重新归一化alpha 取 0.5 这类值而不是随意加一个常数。5.4 带宽选择让 KDE 方向性消失现象KDE 版本算出来两个方向的 TE 几乎一样或者数值大到离谱。原因bandwidth 设得太小时KDE 密度函数退化成一系列脉冲每个样本点只看得到自己条件概率比趋近于 1TE 被高估bandwidth 设得太大时所有密度函数都变成同一个扁平高斯条件概率比趋近于 1TE 又被低估。两个极端都会抹平方向差异。解决用 Scott 规则给一个起点bandwidth N ** (-1 / (d 4))其中 d 是维度。三维联合密度用 d3二维边际用 d2但 KDE 函数里我用的是同一个带宽简化处理。靠谱的做法是分别用 0.5、1、2 倍默认带宽各跑一遍看 TE 峰值位置是否稳定峰值位置能稳定数值有差别结论才值得信任。5.5 不做置换检验任何正数 TE 都没有意义现象算出的 TE 0.02看起来是正数于是宣布找到了方向性耦合。原因传递熵的估计量是正偏的直方图法和 KDE 法在有限样本下都会产生正的“基础噪声”。两个白噪声序列也能算出一个明显大于零的 TE这并不代表有信息流。解决做零假设检验。把发送方序列循环平移随机步破坏原始配对关系再重新计算 TE重复 200 到 1000 次得到零分布。如果真实 TE 超过了零分布的 95% 分位点才可以说存在显著信息流。循环平移比简单打乱好因为它能保留发送方序列内部的时间相关性避免因为自相关结构改变而产生假结论。6. 进阶用法多变量条件 TE 与稳定性验证6.1 条件传递熵剔除公共驱动当第三个序列 Z 同时驱动 X 和 Y 时普通 TE 会把虚假方向算出来。条件传递熵在估计条件概率时把 Z 的过去也加入条件集合也就是把p(x_{ns} | x_n, y_n)改成p(x_{ns} | x_n, y_n, z_n)并与不包含 y_n 的条件概率做比值。实现上就是在直方图或 KDE 里多加一列数据样本需求会大幅上升所以高阶场景建议直接用符号化。6.2 符号化排序近似符号化 TE 不估计连续密度而是把三维样本按数值排序。对每个时刻 t看x_past、y_past、x_fut的大小顺序用排列模式代替原始值然后在离散模式上计数计算条件互信息。这个方法对异常值和噪声极稳参数只有一个延迟 s不需要调带宽。缺点是三变量的排序模式最多 6 种信息分辨率有限适合作为 TE 结果的交叉验证。6.3 稳定性验证我自己的三条结局检查我现在拿到任何 TE 结果都要求自己把这三件事做齐第一画 TE 随延迟变化的曲线而不是只报一个峰值第二把数据切成两段分别计算看峰值延迟是否一致第三至少做 500 次置换检验把原始 TE 和零分布画在同一张图上。如果峰值移动超过两个采样点或者 p 值大于 0.05我宁可不出报告。多年下来这个习惯帮我抓住了很多本来会写进结论里的假因果。实现层面还有一个实用技巧把延迟扫描和带宽测试封装成一个循环输出一个三列数组包含延迟、两个方向的 TE、以及置换 p 值这样每个项目只要跑一次结果就能直接进周报和幻灯片。做传递熵计算时“算得出来”永远不是终点“稳定可复现”才是。希望帮到你。本文还有配套的精品资源点击获取