别再怕‘熵’了！从信息论到脑网络：用Python可视化讲透传递熵与相位传递熵

用Python揭开脑网络信息流动的奥秘传递熵与相位传递熵实战指南神经科学研究中理解不同脑区之间的信息流动一直是核心挑战。传统方法往往依赖复杂的数学公式让许多研究者望而生畏。本文将用Python可视化工具带您直观理解信息论中的熵概念并实战演示如何用传递熵和相位传递熵分析脑网络信息流向。1. 从日常生活理解熵的概念熵在物理学和信息论中有着不同的含义但核心思想都是描述不确定性或混乱程度。想象一下您的办公桌当所有文件整齐归类时您能快速找到需要的文档低熵状态而当文件杂乱堆积时寻找特定文档就变得困难高熵状态。在信息论中香农熵量化了一个随机变量的不确定性。公式虽然重要但更重要的是理解其直观意义import numpy as np def shannon_entropy(prob_dist): 计算香农熵 return -np.sum(prob_dist * np.log2(prob_dist)) # 示例硬币投掷的熵 fair_coin np.array([0.5, 0.5]) # 公平硬币 biased_coin np.array([0.9, 0.1]) # 有偏见的硬币 print(f公平硬币的熵: {shannon_entropy(fair_coin):.2f} bits) print(f偏见硬币的熵: {shannon_entropy(biased_coin):.2f} bits)执行结果公平硬币的熵: 1.00 bits 偏见硬币的熵: 0.47 bits这个简单的例子展示了结果越不可预测公平硬币熵值越高结果越可预测偏见硬币熵值越低。2. 信息论中的进阶概念从联合熵到传递熵理解单个变量的熵后我们需要扩展到多变量系统。以下是关键概念的直观解释联合熵描述两个变量组合的不确定性条件熵已知一个变量后另一个变量的剩余不确定性互信息两个变量共享的信息量这些概念的关系可以用Python可视化import matplotlib.pyplot as plt from matplotlib_venn import venn2 plt.figure(figsize(8, 4)) # 韦恩图展示信息关系 plt.subplot(121) venn2(subsets(0.3, 0.2, 0.1), set_labels(X, Y)) plt.title(互信息可视化) # 信息流动示意图 plt.subplot(122) plt.arrow(0.2, 0.5, 0.6, 0, head_width0.05, head_length0.05, fcr) plt.text(0.5, 0.6, 信息流向, hacenter) plt.xlim(0, 1) plt.ylim(0, 1) plt.axis(off) plt.title(信息流动方向) plt.tight_layout() plt.show()传递熵在此基础上更进一步它能够识别信息流动的方向性。在脑网络分析中这特别有价值——我们不仅能知道两个脑区是否相关还能知道信息主要从哪个脑区流向哪个脑区。3. 传递熵的Python实现与可视化传递熵的核心思想是如果信号X的历史能帮助预测信号Y的未来超过Y自身历史所做的预测那么存在从X到Y的信息流。以下是计算传递熵的关键步骤将连续信号离散化使用直方图方法计算各种联合概率分布根据传递熵公式组合这些概率from scipy.stats import entropy from sklearn.preprocessing import KBinsDiscretizer def transfer_entropy(x, y, delay1, n_bins5): 计算从x到y的传递熵 # 离散化时间序列 discretizer KBinsDiscretizer(n_binsn_bins, encodeordinal, strategyuniform) x_d discretizer.fit_transform(x.reshape(-1, 1)).flatten().astype(int) y_d discretizer.fit_transform(y.reshape(-1, 1)).flatten().astype(int) # 创建历史序列 y_past y_d[delay:-delay] y_future y_d[delay1:] x_past x_d[:-2*delay] # 计算各种概率分布 p_yf_yp, _ np.histogramdd(np.vstack([y_future, y_past]).T, binsn_bins) p_yp, _ np.histogram(y_past, binsn_bins) p_yp_xp, _ np.histogramdd(np.vstack([y_past, x_past]).T, binsn_bins) p_yf_yp_xp, _ np.histogramdd(np.vstack([y_future, y_past, x_past]).T, binsn_bins) # 归一化为概率 p_yf_yp p_yf_yp / np.sum(p_yf_yp) p_yp p_yp / np.sum(p_yp) p_yp_xp p_yp_xp / np.sum(p_yp_xp) p_yf_yp_xp p_yf_yp_xp / np.sum(p_yf_yp_xp) # 计算各项熵值 H_yf_yp entropy(p_yf_yp.flatten()) H_yp entropy(p_yp) H_yp_xp entropy(p_yp_xp.flatten()) H_yf_yp_xp entropy(p_yf_yp_xp.flatten()) # 传递熵 TE H_yf_yp H_yp_xp - H_yp - H_yf_yp_xp return TE # 示例使用 np.random.seed(42) x np.random.normal(size1000) y np.roll(x, 1) 0.1*np.random.normal(size1000) # y滞后于x te_xy transfer_entropy(x, y) te_yx transfer_entropy(y, x) print(fX→Y的传递熵: {te_xy:.4f}) print(fY→X的传递熵: {te_yx:.4f})执行结果可能类似于X→Y的传递熵: 0.1523 Y→X的传递熵: 0.0021这个结果清楚地显示了信息主要从X流向Y这与我们构造数据的方式一致y依赖于x的历史值。4. 相位传递熵分析脑网络信息流的利器相位传递熵是传递熵在信号相位上的应用特别适合分析神经振荡信号如EEG、MEG。它的计算步骤包括通过希尔伯特变换提取信号相位将相位值离散化计算传递熵from scipy.signal import hilbert def phase_transfer_entropy(x, y, delay1, n_bins8): 计算相位传递熵 # 希尔伯特变换提取相位 phase_x np.angle(hilbert(x)) phase_y np.angle(hilbert(y)) # 将相位从[-π,π]转换到[0,2π] phase_x (phase_x np.pi) % (2*np.pi) phase_y (phase_y np.pi) % (2*np.pi) # 离散化相位 bins np.linspace(0, 2*np.pi, n_bins1) x_d np.digitize(phase_x, bins) - 1 y_d np.digitize(phase_y, bins) - 1 # 计算传递熵复用前面的函数 return transfer_entropy(x_d.reshape(-1,1), y_d.reshape(-1,1), delay, n_bins) # 生成示例信号 t np.linspace(0, 10, 1000) x np.sin(2*np.pi*5*t) 0.5*np.random.normal(size1000) y 0.8*np.sin(2*np.pi*5*(t-0.02)) 0.5*np.random.normal(size1000) # y滞后于x pte_xy phase_transfer_entropy(x, y) pte_yx phase_transfer_entropy(y, x) print(fX→Y的相位传递熵: {pte_xy:.4f}) print(fY→X的相位传递熵: {pte_yx:.4f}) # 计算标准化PTE (dPTE) dPTE_xy pte_xy / (pte_xy pte_yx) print(fX→Y的标准化相位传递熵: {dPTE_xy:.4f})典型输出可能为X→Y的相位传递熵: 0.3214 Y→X的相位传递熵: 0.1025 X→Y的标准化相位传递熵: 0.7582标准化相位传递熵(dPTE)的值在0到1之间dPTE 0.5 表示信息主要从X流向YdPTE 0.5 表示信息主要从Y流向XdPTE ≈ 0.5 表示双向信息流动平衡5. 脑网络信息流可视化实战有了传递熵和相位传递熵的计算方法我们可以将其应用于多通道脑电数据构建完整的信息流网络。以下是关键步骤预处理EEG数据滤波、去噪计算所有通道对的dPTE值构建有向加权网络可视化信息流网络import numpy as np import matplotlib.pyplot as plt import networkx as nx # 模拟5通道EEG数据 (1000时间点) np.random.seed(42) n_channels 5 n_samples 1000 data np.random.normal(size(n_samples, n_channels)) # 添加一些耦合关系通道1→通道2通道3→通道4 data[:,1] 0.5*np.roll(data[:,0], 2) # 通道1影响通道2延迟2个样本 data[:,3] 0.3*np.roll(data[:,2], 3) # 通道3影响通道4延迟3个样本 # 计算所有通道对的dPTE矩阵 dPTE_matrix np.zeros((n_channels, n_channels)) for i in range(n_channels): for j in range(n_channels): if i ! j: pte_ij phase_transfer_entropy(data[:,i], data[:,j]) pte_ji phase_transfer_entropy(data[:,j], data[:,i]) dPTE_matrix[i,j] pte_ij / (pte_ij pte_ji 1e-10) # 避免除以零 # 创建有向图 G nx.DiGraph() nodes [fCh{i1} for i in range(n_channels)] G.add_nodes_from(nodes) # 添加边只保留显著的信息流 threshold 0.6 # 只显示dPTE0.6的连接 for i in range(n_channels): for j in range(n_channels): if i ! j and dPTE_matrix[i,j] threshold: G.add_edge(nodes[i], nodes[j], weightdPTE_matrix[i,j]) # 可视化 plt.figure(figsize(10, 8)) pos nx.circular_layout(G) edge_weights [G[u][v][weight]*3 for u,v in G.edges()] # 线宽反映dPTE大小 nx.draw_networkx_nodes(G, pos, node_size800, node_colorlightblue) nx.draw_networkx_edges(G, pos, widthedge_weights, edge_colorred, arrowsize20) nx.draw_networkx_labels(G, pos, font_size12, font_weightbold) # 添加边权重标签 edge_labels {(u,v): f{G[u][v][weight]:.2f} for u,v in G.edges()} nx.draw_networkx_edge_labels(G, pos, edge_labelsedge_labels) plt.title(脑网络信息流方向 (dPTE), fontsize14) plt.axis(off) plt.tight_layout() plt.show()这张可视化图会清晰地显示哪些脑区之间存在显著的信息流动以及流动的方向和强度。在实际研究中这种分析可以帮助识别脑网络中的信息枢纽和关键传导路径。