简介这是一份用于计算两个时间序列间双向传递熵Transfer Entropy的MATLAB脚本资源面向数据科学、复杂系统分析及神经科学、金融等领域研究人员帮助量化信息流动方向与强度。脚本基于信息熵理论通过预处理时间序列、估计概率分布、计算条件概率等步骤输出A→B与B→A两个方向的传递熵值同时支持在widelymfx框架下扩展应用便于揭示系统内部隐藏的动态交互规律。压缩包为RAR格式仅含1个m文件大小1KB属于精简的独立功能脚本阅读和复用成本低。当前已有453人学习适合需要快速实现传递熵计算、验证算法原理或嵌入自建分析流程的中高级用户选用。1. 传递熵transfer entropy相关不等于因果它补上的是“方向”这块拼图两个变量相关系数 0.85你能说清是谁带动谁吗回归系数再显著也分不清是 X 带 Y 还是 Y 带 X——如果背后还藏着一个共同驱动源结论说反几乎是常态。transfer_entropy传递熵就是冲这个来的它给每对时间序列算一个有方向的数值告诉你“已知目标历史之后额外知道源变量的历史能让未来预测的不确定性降低多少”。它不要求线性假设也不受对称性限制适合金融资产溢出效应分析、脑电肌电耦合、工业传感器故障传播以及建模前的因果预筛。对从业者来说它的价值不是替代回归而是让“先问方向、再论强度”这件事有一个非参数的工具。2. 先看懂传递熵的计算逻辑从熵到条件互信息再到有方向的度量传递熵听起来抽象但拆开看只有三层熵、条件熵、方向。先不急着贴公式我们从信息量本身推过去。2.1 熵、条件熵和互信息信息论里的“不确定性”记账本信息熵衡量一个随机变量的不确定性H(X) -Σ p(x) log2 p(x)单位是 bit。熵越大越难预测。如果我们已经知道了 YX 还剩多少不确定性就叫条件熵 H(X|Y)。互信息 MI(X;Y) H(X) - H(X|Y)说的是“知道 Y 之后 X 的不确定性降低了多少”。互信息有两个特点它总是大于等于 0它是对称的MI(X;Y) MI(Y;X)。这带来一个工程上的尴尬——我们做时间序列分析时最想知道的恰恰是方向到底是 A 影响 B还是 B 影响 A。互信息给不出这个答案。而且它没有时间结构你把 X 和 Y 所有时刻的观测混在一起算它只能告诉你整体关联不能告诉你“A 昨天的值对 B 今天的值有没有额外解释力”。2.2 传递熵的定义在条件里多加一个“源变量的过去”传递熵是在条件互信息框架上加了时间轴。以最常用的一阶形式为例计算 S → X 的传递熵TE(S→X) Σ p(x_{t1}, x_t, s_t) · log2 [ p(x_{t1} | x_t, s_t) / p(x_{t1} | x_t) ]逐个符号解释。x_t 是目标变量 X 在 t 时刻的值x_{t1} 是下一时刻的值s_t 是源变量 S 在 t 时刻的值。整个式子读出来就是在已经知道 X 自己历史 x_t 的情况下如果额外加上 S 的当前值 s_t对 x_{t1} 的预测不确定性还能降低多少。如果 S 对 X 没有真正的影响分子和分母几乎一样TE 接近 0如果有影响这个比值大于 1取对数后得到一个正的 bit 数。严格来说这个式子等价于条件互信息 I(x_{t1}; s_t | x_t)前提是我们只使用一阶历史。如果你要把历史拉长到 k 步和 l 步公式就变成TE(S→X) Σ p(x_{t1}, x_t^{(k)}, s_t^{(l)}) · log2 [ p(x_{t1} | x_t^{(k)}, s_t^{(l)}) / p(x_{t1} | x_t^{(k)}) ]其中 x_t^{(k)} 表示 X 过去 k 个时刻的联合取值s_t^{(l)} 同理。这个更一般的形式在实际代码里很少直接实现因为维度会爆炸后面第 4 章会讲怎么绕开。注意一个关键性质TE(S→X) 和 TE(X→S) 不相等。这是因为条件项里各自用的是自己的历史分子分母拆开之后方向不对称。这正是它区别于互信息的地方也是你判断“谁影响谁”的依据。2.3 和 Granger 因果、相关系数的边界什么时候用传递熵从业者最常问的是我直接用 Granger 因果不行吗线性回归跑一下F 检验一上不就有方向了答案是线性场景下 Granger 因果很好用但有两个硬边界。第一Granger 因果本质是线性自回归模型的残差比较。它检验的是“加入 S 的历史后X 预测残差的方差是否显著下降”。如果真实耦合是非线性的比如阈值效应、符号依赖、状态切换线性模型可能完全测不出来或者测出来的方向是错的。传递熵不假设任何函数形式它直接对比概率分布非线性关系照样进到联合分布里。第二Granger 因果对平稳性、方差齐性有较强依赖。金融收益率序列的波动聚集、工业传感器的工况切换都会破坏这些假设。传递熵的经验分布估计对非平稳更敏感但它的敏感方式不同——后面第 5 章会专门讲非平稳的处理而不是说它天生免疫。至于相关系数它连方向都没有而且只能捕捉线性共变关系。我的习惯是把三者放在一起用先用相关系数做初步筛选再用 Granger 因果看线性层面的方向最后用传递熵补非线性部分。三者一致时结论才敢写进报告不一致时优先怀疑数据里有共同驱动源或滞后阶数没对齐。3. 用 Python 零依赖算传递熵分箱、滞后与一个最小可复现实现理论立住之后落地只需要 NumPy。下面这个实现不依赖任何专业信息论库所有逻辑透明可改适合你拿自己的数据跑第一版结果。3.1 分箱把连续信号变成离散符号传递熵的概率估计需要一个离散空间。常见做法有两种等宽分箱和分位数分箱。等宽分箱按数值范围均匀切但金融数据和传感器数据都有重尾等宽分箱经常出现“一个箱子里装了 90% 的点其它箱子全是空的”这种局面。我一般用分位数分箱把数据按排序后的分位点切保证每个箱子里样本量接近。这样概率估计的稳定性好很多代价是极端值的具体数值被抹掉了但我们算的是信息流动不是数值预测这个代价可以接受。3.2 最小实现代码用频率估计联合概率下面是一份可以直接保存使用的实现核心逻辑不到五十行import numpy as np def _quant_bin(x: np.ndarray, n_bins: int) - np.ndarray: 分位数分箱返回 0..n_bins-1 的整数符号序列。 edges np.quantile(x, np.linspace(0.0, 1.0, n_bins 1)) edges[0] - 1e-9 # 保证最小值落入第一个箱 edges[-1] 1e-9 # 保证最大值落入最后一个箱 return np.digitize(x, edges[1:-1]) def te(src: np.ndarray, tgt: np.ndarray, lag: int 1, n_bins: int 5) - float: 计算 src - tgt 的传递熵单位 bit。 一阶马尔可夫假设预测 tgt[t1] 时只用 tgt[t] 和 src[t]。 s _quant_bin(src, n_bins) x _quant_bin(tgt, n_bins) n len(x) if n lag 10: raise ValueError(序列太短联合分布估计不可靠) # 对齐s[t], x[t] 预测 x[t1] s_past s[lag - 1:] # 源在相对时刻 t 的取值 x_past x[lag - 1:-1] # 目标在 t 的取值 x_future x[lag:] # 目标在 t1 的取值 # 频次统计这是频率估计概率的全部基础 joint np.zeros((n_bins, n_bins, n_bins)) # (s, x_past, x_future) cond_sx np.zeros((n_bins, n_bins)) # (s, x_past) cond_x np.zeros((n_bins, n_bins)) # (x_past, x_future) marg_x np.zeros(n_bins) # (x_past) for i in range(len(x_past)): joint[s_past[i], x_past[i], x_future[i]] 1 cond_sx[s_past[i], x_past[i]] 1 cond_x[x_past[i], x_future[i]] 1 marg_x[x_past[i]] 1 n_eff len(x_past) out 0.0 eps 1e-9 for si in range(n_bins): for xi in range(n_bins): for xj in range(n_bins): cnt joint[si, xi, xj] if cnt 0: continue p_joint cnt / n_eff # p(x_future | s, x_past) p_full cnt / (cond_sx[si, xi] eps) # p(x_future | x_past) p_base cond_x[xi, xj] / (marg_x[xi] eps) out p_joint * np.log2(p_full / p_base eps) return out代码拆开看就三件事分箱、数频次、按公式求和。joint 数组三个维度分别对应源当前值、目标当前值、目标未来值cond_sx 统计的是在给定源和目标当前值的条件下有多少样本cond_x 和 marg_x 用于计算基线预测概率。log 以 2 为底输出单位是 bit这样不同实验之间可以横向比较。lag 参数在这里的语义要特别说清楚它表示源序列和目标序列之间的时间偏移。代码里通过 s_past s[lag - 1:] 实现——lag1 时用的是当前时刻的源值lag2 时用的是上一时刻的源值。如果你确信业务机制是“昨天的 S 影响今天的 X”就把 lag 设为对应步数而不是机械地写 1。n_bins 我建议从 5 开始数据量低于一千点时不要超过 6后面第 4 章会展开讲。3.3 合成数据验证先确认它能恢复正确的方向拿到代码第一件事不是跑真实数据而是造一份已知耦合方向的合成序列确认方向检测不是玄学。下面这个例子构造了一个带自回归和外部输入的序列rng np.random.default_rng(7) T 5000 src rng.standard_normal(T) # 源纯白噪声 noise rng.standard_normal(T) tgt np.zeros(T) for i in range(2, T): # S 延迟 2 期影响 T同时 T 有自相关 tgt[i] 0.55 * tgt[i - 1] 0.6 * src[i - 2] 0.3 * noise[i] print(fsrc - tgt: {te(src, tgt, lag2, n_bins6):.4f} bit) print(ftgt - src: {te(tgt, src, lag2, n_bins6):.4f} bit)跑出来的结果里src - tgt 应该明显大于 tgt - src。如果两个方向都接近 0先检查数据量再做分箱和滞后检查如果反向大于正向多半是共同驱动或非平稳导致的假象对应第 5 章的避坑条目。每次换数据集我都走一遍这个合成验证流程确保代码链路本身没毛病再去看业务结论。4. bins、lag 和显著性怎么定三个必调参数与一套验证方法传递熵最让人头疼的不是公式而是三个参数的设置分箱数、滞后阶数、以及算出来的数值到底算不算显著。这三个问题不解决TE 就只是个黑匣子。4.1 bins 选多少五分位起步扫一遍看稳定性bins 太少信息被过度离散化真实的信息流动可能被抹平bins 太多联合分布有 n_bins^3 个格子数据量稍微不够就全是零频频率估计的方差急剧变大。实践中我以 5 为起点数据量每翻十倍可以加一档最多用到 8。更可靠的做法是跑一个 bins 扫描把结论稳定性直接摆出来看for b in [4, 5, 6, 8, 10]: te_forward te(src, tgt, lag2, n_binsb) te_backward te(tgt, src, lag2, n_binsb) print(fbins{b:2d} forward{te_forward:.4f} bit backward{te_backward:.4f} bit)如果不同 bins 下方向结论始终一致这个结论才是可信的。如果 bins4 说正向bins8 又说反向那大概率是数据量不足或存在共同驱动而不是 bins 本身的问题。我在项目里遇到这种情况会直接停掉 TE 分析先去查数据生成过程。4.2 lag 怎么找按业务周期扫选显著且峰值稳定的滞后lag 设错是新手最容易翻车的地方。真实耦合延迟两期你却用 lag1算出来的根本不是延迟因果而是同一时刻的同步互信息方向会变得没有意义。我的做法是先按业务经验缩小范围日频金融数据试 1 到 5 天传感器数据按采样周期试对应的时间窗口然后做 lag 扫描for lag in range(1, 7): te_fwd te(src, tgt, laglag, n_bins6) te_bwd te(tgt, src, laglag, n_bins6) print(flag{lag}: forward{te_fwd:.4f} bit backward{te_bwd:.4f} bit)正确设置下TE 曲线通常会在真实延迟处出现峰值之后缓慢衰减。如果所有 lag 都差不多大说明可能没有方向性耦合也可能是共同驱动源在起作用。有一点要强调不要只看 TE 最大值的 lag 就下结论要配合下一节的置换检验确认峰值是统计显著的。4.3 置换检验算出 0.05 bit 到底是信号还是噪声TE 是一个非负统计量纯粹由噪声驱动的序列也能算出正数。所以必须回答一个问题这个正数有没有超过随机水平标准做法是有效传递熵和置换检验把源序列在时间上打乱重新计算 TE重复多次得到一个零假设分布。def te_surrogates(src, tgt, lag1, n_bins5, n_surr200, seed42): rng np.random.default_rng(seed) obs te(src, tgt, lag, n_bins) surr [] for _ in range(n_surr): src_shuf rng.permutation(src) # 打乱源序列顺序 surr.append(te(src_shuf, tgt, lag, n_bins)) p_value (1 sum(1 for s in surr if s obs)) / (1 n_surr) return obs, np.percentile(surr, 95), p_value需要提示一个细节直接整体打乱源序列会破坏它自身的自相关结构导致置换分布失真。严格做法是块置换把序列切成若干连续块再整块重排保序内部自相关。数据量在几千点以上时我用块长约为 lag 的十倍数据量小就放宽到直接置换并注明代价。输出里我会同时记录原始 TE 和 95% 分位数如果原始值低于百分位就直接标记为不显著不再用来做业务判断。这个流程在多变量场景也同样适用。标题里的 widelymfx 对应的就是金融多资产交叉分析把每对变量的 TE 都算出来组织成交叉矩阵再做置换检验过滤掉噪声边。技术栈上既有 R 生态封装好的现成包也有人用 Python 自己拼矩阵核心都是先过显著性这一关矩阵才有意义。5. 传递熵避坑指南五条让我翻过车的高频问题这一章全部来自真实踩坑记录。每一条都按“现象 → 原因 → 解决”写照着排查比自己瞎调快得多。5.1 现象单向耦合算成双向方向结论全反合成数据明明只有 S 影响 XTE 却给出两个方向都显著甚至反向更大。原因最常见的凶手是共同驱动源。S 和 X 都受一个未观测变量 Z 影响时即便 S 对 X 没有任何直接作用TE(S→X) 也会显著因为 S 携带了 Z 的信息这部分信息对预测 X 有贡献。其次是数据非平稳两个序列有相同的趋势项趋势本身就会制造虚假的预测关系。解决先做平稳化处理差分或去趋势后再算如果业务场景里确实存在共同驱动变量改用第 6 章的条件传递熵把它放进条件项控制住。判断共同驱动的一个快速方法是看两个序列的互信息是否远大于两个方向的 TE 最大值——如果是优先怀疑 Z。5.2 现象换个 bins 值结论就翻转bins5 时正向显著bins8 时反向显著bins3 时啥都不显著。原因数据量撑不起高维联合分布。n_bins8 时联合空间有 512 个格子一千点数据平均一个格子只有两个样本频率估计全是噪声这时候左右方向的 TE 差没有意义。反过来 bins 太小又把动态模式压缩掉了。解决先做第 4.1 节的 bins 扫描强制看结论稳定性。数据量少于五千点bins 不要超过 6如果换 bins 后方向翻转把这次分析视为数据不足而不是参数没调好。多攒数据或者降采样都行就是别试图靠调 bins 调出想要的结论。5.3 现象lag 设成 1把同步关系当成了因果两个变量其实是同向联动没有谁导谁但 lag1 时 TE 大得吓人。原因当两个序列的当期值高度同步时用 s[t] 预测 x[t1] 很可能比只用 x[t] 预测更好因为 s[t] 实际上预演了 x[t]当前已发生但还没进模型的信息。这是典型的“未来信息泄漏”——不是真正的传递。解决把 lag 从 1 逐步加到 5 以上观察 TE 曲线。如果是同步假象TE 在 lag1 时高lag 加大后迅速跌到接近零如果是真实延迟因果TE 会在匹配真实延迟的地方出现峰值。另外算之前把两个序列做滞后交叉相关性图也有助于初判真实延迟范围。5.4 现象非平稳序列跑 TE数值忽大忽小同一对变量换一段时间的窗口TE 从 0.2 bit 变成 0.02 bit方向也变。原因非平稳导致边缘分布漂移分位数分箱的边界随窗口变化同样的耦合强度在不同窗口里被编码成不同的符号序列TE 数值自然不稳定。如果序列有突变或断点频率估计更是直接失真。解决先做差分或提取残差确保进入 TE 的是平稳序列。金融数据用对数收益率而不是价格传感器数据先去除工况切换造成的基线漂移。对长周期数据我习惯用滚动窗口分别计算 TE见第 6.2 节观察数值是否随时间漂移漂移本身也是一种结论但不要把漂移的原始数据整合成一个 TE 值写进报告。5.5 现象用 TE 绝对值大小直接排名和业务常识对不上TE 排序说 A 对 B 的影响最强业务经验里明明是 C 主导。原因TE 的绝对数值受 bins、lag、数据量、噪声水平影响不同变量对的“底噪”不一样数值小不代表没影响数值大也可能是噪声抬高。把 TE 当强度比大小是不对的它能给的是方向和显著性不是标准化的效应量。解决横向比较时统一所有参数配置包括相同的 bins、lag、数据长度和平稳化流程然后只比较“显著边”的相对排名不显著的边一律不参与讨论最后一定用置换检验得到 p 值而不是拿原始 TE 排序。报告里我会同时列原始 TE 和 p 值让读者自己判断强弱。6. 多变量交叉场景的动态信息网络从条件传递熵到滚动验证单变量 TE 只是起步。真实场景里变量之间互相牵制条件传递熵和滚动窗口才是能直接落到业务上的进阶手段。6.1 条件传递熵把共同驱动源控制住条件传递熵在条件项里除了目标变量自身历史再加入一个或多个控制变量 Z公式写成CTE(S→X | Z) Σ p(x_{t1}, s_t, x_t, z_t) · log2 [ p(x_{t1} | s_t, x_t, z_t) / p(x_{t1} | x_t, z_t) ]它回答的问题是在已经知道 Z 的前提下S 对 X 还有没有额外的信息贡献。实现上只需要在频次统计里多加一个维度但注意维度从三维变成四维数据量需求大幅上升数据不够时结果会比单变量 TE 更不稳定。我的处理顺序是先跑单变量 TE 做初筛找到可能的因果关系候选再用条件 TE 验证关键边而不是一上来就对所有组合跑条件版本。6.2 滚动窗口与交叉矩阵一份可照抄的最小工作流拆解一个完整可复用的流程对每个窗口内的数据逐对计算 TE 和置换检验 p 值组成方向矩阵然后观察矩阵随时间的变化def rolling_te_matrix(data: list, window: int, step: int, lag: int, n_bins: int): n_var len(data) results [] for start in range(0, len(data[0]) - window, step): mat np.zeros((n_var, n_var)) for i in range(n_var): for j in range(n_var): if i j: continue obs, th95, p te_surrogates( data[i][start:start window], data[j][start:start window], laglag, n_binsn_bins, n_surr100) mat[i, j] obs if p 0.05 else 0.0 results.append(mat) return results这份代码把“不显著就置零”作为默认策略是我所有 TE 项目的基本盘。业务上可以进一步做聚类、找网络中枢但核心方向判断在前一步就已经完成。后来我所有 TE 分析都强制过一遍这个流程——合成验证、参数扫描、置换检验、滚动观察——基本没有再翻过车。希望帮到你。本文还有配套的精品资源点击获取