传递熵原理与MATLAB实现：从信息论到延迟时间因果分析

简介面向信息论与复杂系统分析研究者的MATLAB实现围绕传递熵计算与最优延迟时间估计展开。程序通过联合概率分布与条件熵计算量化两个系统之间信息传递强度并扫描候选延迟以定位使传递熵最大的时间点为因果推断和短期预测提供可复现的数值工具。压缩包内共3个文件包含2个.m脚本与1张结果示意图整体仅20KB结构精简适合已有概率统计基础的研究生或工程师快速阅读与扩展。目前已有221人学习示意图可作为运行结果的对照参考。除数值实现外资源还体现了最大熵原理在缺乏先验约束时如何选取最不确定分布这一核心思想理解延迟时间对传递熵峰值的影响后可将该方法迁移至通信网络效能评估、生物信号耦合分析和金融时序关联等研究场景。1. 传递熵计算先搞清楚它在算什么如果手上有两条时间序列比如脑电不同通道的信号或者金融里两个板块的指数你想回答的通常不是「它们有没有关系」而是「谁在驱动谁滞后多久」。相关系数和互信息只能给出对称的关联度分不清方向传递熵Transfer Entropy则专门补这个缺口——它衡量的是在已知目标序列历史信息的前提下额外知道源序列的历史能让目标未来状态的不确定性降低多少。压缩包里的传递熵传递时间计算.m就是一个用 MATLAB 实现传递熵估计的程序并且它把「延迟时间」也纳入了计算流程通过扫描不同的时间偏移找到传递熵最大的那个迟滞点。这个资源适合做时间序列因果分析的人无论是信号处理、金融数据分析还是复杂系统的实验数据处理拿到手后可以直接改数据路径跑起来关键是搞清楚它内部怎么算后面调参数才不会翻车。2. 传递熵与最大熵公式拆解和概率估计的底层逻辑2.1 从熵到传递熵一条公式理清信息流方向信息熵衡量不确定性公式是 H(X) -Σ p(x) log p(x)。如果知道了另一变量 Y条件熵 H(X|Y) -Σ p(x,y) log p(x|y) 表示剩余不确定性H(X) 减去 H(X|Y) 得到的互信息刻画共享信息但这个值是对称的——它不告诉你信息是 X 流向 Y还是 Y 流向 X。传递熵打破对称性。把两个时间序列记为 X_t 和 Y_t我们考虑从 X 到 Y 的有向信息传递。设 Y 的未来是 y_{t1}Y 的历史是 y_t^{(k)}X 的历史是 x_t^{(l)}则 TE_{X→Y} 定义为TE_{X→Y} H(y_{t1} | y_t^{(k)}) - H(y_{t1} | y_t^{(k)}, x_t^{(l)})展开成联合概率形式就是TE_{X→Y} Σ p(y_{t1}, y_t^{(k)}, x_t^{(l)}) log [ p(y_{t1} | y_t^{(k)}, x_t^{(l)}) / p(y_{t1} | y_t^{(k)}) ]这个式子的含义很直接如果加入 X 的历史之后Y 未来状态的条件熵明显下降说明 X 对 Y 存在可检测的信息传递。反过来把 X 和 Y 对调得到 TE_{Y→X}比较两个方向的传递熵就能判断主导关系。.m文件里后半部分所做的工作本质上就是在估计这个比值里的各个概率项。下面这张表列出了程序里会出现的主要符号后面读代码时可以直接对照符号含义典型取值X, Y源序列和目标序列列向量长度 N≥200tau延迟步数1~50按采样率定bins直方图分箱数4~8alpha拉普拉斯平滑系数0.01 或 1H1条件熵 H(y_future | y_past)正数单位 nat 或 bitH2条件熵 H(y_future | y_past, x_past)小于等于 H12.2 最大熵原理不知道概率时怎么猜才不偏传递熵估计最大的难点是联合概率 p(y_{t1}, y_t^{(k)}, x_t^{(l)})。样本有限时你永远不知道真实分布长什么样。最大熵原理给了一个选择标准在没有额外信息约束的前提下应该选择熵最大的概率分布因为这是最保守、最不会引入人为偏见的假设。具体到实现里最常见的做法是等宽直方图。把每个变量的取值范围切成 bins 个格子统计每个格子的样本频次频次除以总样本数就是联合概率的估计值。这个操作无形中就在贯彻最大熵思想——每个格子的先验概率被设成相等的。样本一大联合空间会变得稀疏很多格子频次为零。直接 log(0) 会得到 NaN所以程序里普遍要做拉普拉斯平滑。下面这段代码展示了联合概率估计的核心逻辑function jointP estimate_joint_pdf(data, bins, alpha) % data: n x 3 矩阵三列分别为 y_future, y_past, x_past % bins: 每维分箱数 % alpha: 平滑系数 n size(data, 1); binIdx zeros(n, 3); for d 1:3 lo min(data(:, d)); hi max(data(:, d)); if hi lo, hi lo eps; end edges linspace(lo, hi, bins 1); binIdx(:, d) discretize(data(:, d), edges, IncludedEdge, right); end jointCounts accumarray(binIdx, 1, [bins bins bins]) alpha; jointP jointCounts / sum(jointCounts(:)); endbins是每维分箱数alpha是平滑常数等价于给每个格子预先放了一点虚拟样本。accumarray把每个样本的离散坐标累加成一个三维计数矩阵最后归一化得到概率。这里 alpha 越大估计出的分布越平坦熵越大越符合最大熵的保守原则alpha 太小则容易过拟合单个噪声样本就可能主导某个格子的概率。我一般会先用 alpha1 跑通流程再对比 alpha0.01 看结果是否稳定。2.3 延迟时间为什么不能拍脑袋信息从 X 传到 Y 不可能瞬时完成。脑电信号跨脑区传导要几十毫秒金融市场价格对消息的反应可能需要几分钟或几个采样周期。因此传递熵公式里的 x_t^{(l)} 不应该是当前时刻的 X而应该是 X 在 t-τ 时刻的值。τ 就是延迟时间程序里也叫 tau。如果 tau 取太小消息还没来得及传到 Y传递熵自然很小取太大X 的历史信息可能已经与当前 Y 无关传递熵会被噪声淹没。更麻烦的是真实系统的延迟往往未知甚至可能随时间变化。所以成熟的做法不是手填一个 tau而是在一定范围内扫描对每个候选 tau 计算一次 TE得到一条 TE(tau) 曲线峰值对应的 tau 就是估计出的延迟。这也是为什么压缩包里那个.m程序会包含两层逻辑——内层是单次传递熵计算外层是延迟扫描。先理解了这一点再打开代码就不容易看晕。3. MATLAB 实现从 .m 文件到延迟时间曲线3.1 文件结构一个 .m 脚本包含哪几个模块压缩包解开后核心文件是传递熵传递时间计算.m另外还有一个示意图无标题.jpg大概率是程序跑出来的 TE-延迟曲线截图用来对照结果。.m文件虽然只有一个但内部通常会按模块分段常见结构是这样模块作用对应变量数据加载读入 Excel / txt / mat 数据rawX, rawY参数定义设置 bins、tau 范围、alphabins, tauList, alpha数据预处理去趋势、标准化、去异常值X, Y传递熵核心函数计算给定 tau 下的 TE 值te transfer_entropy_value(...)延迟扫描循环 tau 得到 TE 数组TE_list, tau_best绘图画 TE-tau 曲线并标注峰值figure, plot, xline我习惯把核心计算单独抽成一个函数而不是揉在脚本里。因为扫描延迟时要调用几十次独立函数方便测试也方便你在里面加 disp 打印中间结果。如果打开发现原文件是脚本式写法你也可以自己改造成函数形式不影响结果。3.2 核心代码直方图估计联合概率与传递熵传递熵核心计算可以拆成两个条件熵。第一项 H(y_future | y_past) 只涉及目标和它自己的历史是二维直方图第二项 H(y_future | y_past, x_past) 是三维联合直方图。代码实现如下function te transfer_entropy_value(X, Y, tau, bins, alpha) % X, Y: 列向量时间序列 % tau: 延迟步数正整数 % bins: 每维分箱数量经验值 4~8 % alpha: 拉普拉斯平滑默认 1 if nargin 5, alpha 1; end N length(X); t (tau 1) : N - 1; % 满足 y(t1), y(t), x(t-tau) 都存在的索引 Y_future Y(t 1); Y_past Y(t); X_past X(t - tau); % 把每列离散成 1..bins 的整数索引 function idx toBins(v) lo min(v); hi max(v); if hi lo, hi lo 1; end edges linspace(lo, hi, bins 1); idx discretize(v, edges, IncludedEdge, right); end iyf toBins(Y_future); iyp toBins(Y_past); ixp toBins(X_past); % H(Y_future | Y_past) C1 accumarray([iyf, iyp], 1, [bins bins]) alpha; P1 C1 / sum(C1(:)); H1 -sum(sum(P1 .* log(P1 ./ sum(P1, 1)))); % H(Y_future | Y_past, X_past) C2 accumarray([iyf, iyp, ixp], 1, [bins bins bins]) alpha; P2 C2 / sum(C2(:)); marginal squeeze(sum(C2, 1)); % 对 y_future 求和得到 p(y_past, x_past) H2 0; for i 1:bins for j 1:bins pz marginal(i, j) / sum(C2(:)); pygz P2(:, i, j) / sum(P2(:, i, j)); H2 H2 - pz * sum(pygz .* log(pygz eps)); end end te H1 - H2; end这段代码分四步走先构造三元组样本再离散化到 bin 索引然后分别统计二维和三维联合频次最后计算两个条件熵的差值。accumarray([iyf, iyp], 1, [bins bins])的作用是建立二维频次表P1 ./ sum(P1, 1)表示在给定 y_past 维度条件下求 y_future 的条件概率。第二项循环里要特别留意eps的引入防止某个条件概率为 0 时 log 计算出 -Inf。这正是最大熵平滑的延续alpha 已经避免了零概率eps 只是再兜一层底。3.3 参数设置bin 数、延迟范围、数据归一化参数对结果影响最大也是最容易出现「玄学结果」的地方。先说 bins每维分箱数直接决定联合空间的体积bins5 时三维联合空间有 125 个格子样本量少于一千就容易出现大量空箱。我通常固定 bins5 或 6然后在 alpha0.01 和 alpha1 之间做敏感性对比。延迟范围 tauList 的选取要看业务背景如果是秒级采样信号真实延迟通常不超过几十个采样点可以取 1 到 N/10如果完全未知就先取 1:50 跑一遍观察峰值是否落在边界上落在边界就扩大范围。数据归一化在直方图法里不是必须的因为等宽分箱对线性缩放不敏感。但如果数据有趋势、突变或异常尖峰分箱边界会被异常值拉得很开导致有效信息挤在少数几个 bin 里。我一般会先做 z-score 标准化再对非平稳序列做一阶差分最后才进入传递熵计算。注意差分会改变序列长度延迟扫描的索引也要对应调整否则会报维度不匹配的错误。4. 延迟时间搜索传递熵峰值怎么找才靠谱4.1 扫描延迟的常规写法确定了核心函数后扫描延迟就是一个循环问题。为了让结果更稳建议把 TE 值存成数组方便绘图和后续显著性检验。基本的扫描代码如下tauList 1:30; TE zeros(size(tauList)); X zscore(rawX); % 预处理后的源序列 Y zscore(rawY); % 预处理后的目标序列 for i 1:length(tauList) TE(i) transfer_entropy_value(X, Y, tauList(i), 5, 1); end [maxTE, idx] max(TE); tau_best tauList(idx); fprintf(最佳延迟: %d, TE %.4f\n, tau_best, maxTE);这个写法时间复杂度是 O(maxTau * bins^3)bins 不大于 6 时10000 个样本跑 30 个延迟大约要几秒到十几秒。如果样本量到百万级三重循环会明显变慢可以先减小 bins 到 4或者对原始序列降采样。还有一个实用的提速技巧transfer_entropy_value里的分箱边界只取决于数据范围与 tau 无关。所以可以提前算好所有变量的全局分箱边界再在每个 tau 下只做累积计数而不是重新调 linspace。这样能省掉反复排序和边界计算的开销。4.2 峰值判定除了最大值还要看显著性TE(tau) 曲线取最大值是最直观的但真实数据往往有多个局部峰而且可能有噪声导致的假峰。单纯取最大值可能选到随机波动上这就是很多传递熵程序被吐槽「结果不可复现」的主要原因。更可靠的做法是加一个置换检验nPerm 200; TE_perm zeros(nPerm, length(tauList)); for k 1:nPerm X_shuffled X(randperm(length(X))); for i 1:length(tauList) TE_perm(k, i) transfer_entropy_value(X_shuffled, Y, tauList(i), 5, 1); end end threshold quantile(TE_perm, 0.95, 1); significant TE threshold;置换检验的逻辑是把源序列 X 随机打乱破坏 X 与 Y 的真实时序关系然后重新计算 TE。打乱后的 TE 分布代表「没有真实耦合时候的随机水平」。如果某个 tau 下原始 TE 超过了 95% 的置换值才能认为这个峰是显著的。这个检验很吃算力200 次置换 × 30 个 tau可能要跑几分钟。可以先用 50 次置换粗筛确定显著区间后再对区间内的 tau 做细扫描。4.3 可视化输出传递熵—延迟曲线怎么看绘图不仅是给人看更是排查问题的第一道关卡。常规画法如下plot(tauList, TE, b-o, LineWidth, 1.5); hold on; plot(tauList, threshold, r--, LineWidth, 1.2); xlabel(延迟 tau); ylabel(传递熵); legend({TE_{X→Y}, 95% 置换阈值}, Location, best); grid on;拿到这张图先看三件事第一峰值是否远离 tau 范围的边界如果在 30 处还在上升说明真实延迟可能更大得继续扩大范围第二峰值是否明显高于置换阈值如果整条曲线都在阈值以下这组数据里可能根本没有方向性信息传递第三TE 曲线是否光滑如果锯齿特别剧烈多半是 bins 太大或样本量不足需要降 bins 或加数据。我习惯把正向 TE 和反向 TE 画在同一张图用蓝色和红色区分这样谁主导、滞后多少一眼就能看出来。5. 避坑排查传递熵计算常见翻车点与修正方案5.1 四条高频踩坑记录坑 1TE 结果恒等于 0 或接近 0现象不管怎么调延迟TE 值都在 1e-4 以下和没算一样。 原因大概率是 tau 范围没覆盖真实延迟或者源序列和目标序列在给定延迟下根本没有信息交叠。还有一种可能是分箱边界设置不当数据范围被极值拉宽导致大多数样本落在同一个 bin。 解决先做滞后互相关图看哪个 lag 处相关系数绝对值最大再把这个 lag 作为 tau 搜索中心点同时用分位数截断异常值避免极值拉扁直方图。坑 2程序报 NaN 或 Inf现象跑完 TE 数组里出现 NaN或者 log 计算出-Inf。 原因联合频次里某些格子概率为零log(0)产生了未定义值。 解决在transfer_entropy_value里给所有log参数加eps并且确保alpha大于零。如果使用了P2(:, i, j) / sum(P2(:, i, j))分母永远不会是零因为拉普拉斯平滑已经保证所有格子计数至少为 alpha。坑 3bins 稍微一改结论就翻转现象bins4 时 X→Y 显著bins8 时反而 Y→X 显著。 原因样本量不足以支撑高维联合分布。bins8 时三维联合空间有 512 个格子10000 个样本看起来很多但要填满 512 个格子并且让每个格子计数可靠仍然不够。 解决固定 bins5然后把 alpha 作为敏感性参数。最理想的做法是用自助法重复估计 TE看置信区间是否包含 0。如果区间太宽就不要轻易下因果结论。坑 4信号有趋势算出来的传递熵是假的现象两个独立但有缓慢上升趋势的序列TE 却很高。 原因传递熵本身没有内置平稳性检验趋势项被直方图当成了稳定的联合概率结构从而产生虚假的预测力。 解决对数据先做差分或去趋势对残差做传递熵计算。金融数据尤其要注意这一点价格序列必须先转成收益率或对数差分再做 TE 分析。这是我踩过最深的坑没有之一——第一次用在沪深数据上算出两个指数互相驱动后来发现只是共同趋势在作怪。5.2 从中间结果排查数据问题还是程序问题拿到一个不正常的 TE 结果先不要急着改参数而是打印中间量。在transfer_entropy_value里加两个输出H1 和 H2。如果 H1 小于 H2说明程序逻辑有错——加入更多条件信息后熵不可能变大。如果 H1 正常、H2 等于 H1说明 x_past 没有提供任何额外信息数据本身可能没有耦合。可以用下面这段代码快速定位[te, H1, H2] transfer_entropy_value(X, Y, tau, bins, alpha); fprintf(H1 %.4f, H2 %.4f, TE %.4f\n, H1, H2, te); assert(H2 H1 1e-12, 条件熵逻辑异常请检查联合概率统计);当 H2 只比 H1 小一点点时TE 数值会很小此时主要矛盾不是程序 bug而是数据信噪比过低。处理办法是增加样本量、降低噪声或者尝试用符号化传递熵把数据转成 rank 后再分箱它对噪声的鲁棒性比原始值直方图更好。6. 验证技巧用自造耦合数据检验你的传递熵程序6.1 自造数据一个已知耦合的 VAR 模型拿到这份资源后我建议你做的第一件事不是跑真实数据而是用已知答案的数据验证程序。构造一个最简单的单向耦合系统X 是白噪声Y 在延迟 d 步后受 X 影响。这样真实延迟 d 已知如果程序输出峰值不在 d 附近那程序参数一定有问题。生成代码很简单N 10000; d 3; X randn(N, 1); Y zeros(N, 1); for t (d 1):N Y(t) 0.7 * X(t - d) 0.5 * randn; end这里 d3意味着 X 在 t-3 时刻的信息会影响 Y_t。噪声系数 0.5 保证信噪比合理TE 不会小到测不出来也不会大得没有挑战性。6.2 跑通全流程看峰值是否出现在真实延迟用上一章的transfer_entropy_value函数扫描 tau1:10X zscore(X); Y zscore(Y); tauList 1:10; TE_fwd zeros(1, 10); TE_rev zeros(1, 10); for i 1:10 TE_fwd(i) transfer_entropy_value(X, Y, i, 5, 1); TE_rev(i) transfer_entropy_value(Y, X, i, 5, 1); end [~, bestFwd] max(TE_fwd); [~, bestRev] max(TE_rev); fprintf(正向峰值延迟 %d反向峰值延迟 %d\n, bestFwd, bestRev);正确结果是正向 TE 在 tau3 处出现明显峰值反向 TE 应该整体很小且没有稳定峰值。如果这两个条件不满足问题一般出在数据预处理或 bins 选择上。我第一次跑这个验证时反向 TE 在 tau1 处居然比正向还高查了半天发现是 Y 自身的自相关在捣鬼——后来给 Y 去掉自回归项后结果才干净。这提醒我验证程序时生成数据要尽量简单别混入额外耦合。6.3 把这个验证步骤固化到你的测试习惯里从那以后我每次拿到新的传递熵程序或者调整了 bins、alpha 等参数都会先用这段自造数据强制走一遍全流程。生成数据 → 扫描延迟 → 对比正反向 TE → 确认峰值位置整个验证脚本不到两分钟。只有这关过了我才会把程序用到真实数据上否则即使结果看起来合理你也没法分辨它是真实信息传递还是程序 bug。这份传递熵传递时间计算.m的底子不错核心框架完整但参数敏感性偏高希望你用的时候把这套验证步骤一并变成自己的习惯。希望帮到你。本文还有配套的精品资源点击获取