告别理论恐惧:用Python可视化一步步理解‘熵’到‘传递熵’的因果发现逻辑

告别理论恐惧:用Python可视化一步步理解‘熵’到‘传递熵’的因果发现逻辑 告别理论恐惧用Python可视化一步步理解‘熵’到‘传递熵’的因果发现逻辑信息论中的熵概念常常让人望而生畏那些复杂的公式和抽象的定义让许多开发者望而却步。但事实上这些概念完全可以通过直观的可视化和简单的Python代码来理解。本文将带你用Python从最基础的掷骰子例子开始一步步构建对熵、联合熵、条件熵、互信息最终到传递熵的直观理解。1. 从骰子开始理解基础熵概念熵本质上衡量的是一个系统的不确定性。让我们从一个简单的例子开始掷一个公平的六面骰子。import numpy as np import matplotlib.pyplot as plt # 公平骰子的概率分布 fair_dice np.ones(6)/6 # 作弊骰子总是出6 loaded_dice np.array([0, 0, 0, 0, 0, 1]) def entropy(prob_dist): return -np.sum(prob_dist * np.log2(prob_dist, whereprob_dist0)) print(f公平骰子的熵: {entropy(fair_dice):.3f} bits) print(f作弊骰子的熵: {entropy(loaded_dice):.3f} bits)运行这段代码我们会看到公平骰子的熵约为2.585 bits而总是出6的作弊骰子熵为0。这直观地展示了熵的含义公平骰子的结果不确定性最大熵最高而结果确定的作弊骰子熵为零。注意当概率为零时我们约定0×log00这在信息论中是标准做法。我们可以用柱状图更直观地展示不同概率分布的熵plt.figure(figsize(10,4)) plt.subplot(121) plt.bar(range(1,7), fair_dice) plt.title(f公平骰子 (熵{entropy(fair_dice):.3f} bits)) plt.subplot(122) plt.bar(range(1,7), loaded_dice) plt.title(f作弊骰子 (熵{entropy(loaded_dice):.3f} bits)) plt.show()2. 联合熵与条件熵天气与穿衣的例子现在考虑一个更实际的例子天气与穿衣选择的关系。假设我们有以下联合概率分布短袖 (Y0)长袖 (Y1)P(X)晴天 (X0)0.40.10.5雨天 (X1)0.10.40.5P(Y)0.50.51.0我们可以计算各种熵# 定义联合概率分布 joint_prob np.array([[0.4, 0.1], [0.1, 0.4]]) px joint_prob.sum(axis1) # P(X) py joint_prob.sum(axis0) # P(Y) # 计算各种熵 hx entropy(px) # H(X) hy entropy(py) # H(Y) hxy entropy(joint_prob.flatten()) # H(X,Y) # 条件熵 H(Y|X) H(X,Y) - H(X) hy_given_x hxy - hx print(fH(X): {hx:.3f} bits) print(fH(Y): {hy:.3f} bits) print(fH(X,Y): {hxy:.3f} bits) print(fH(Y|X): {hy_given_x:.3f} bits)这个例子展示了条件熵的含义知道天气情况后穿衣选择的不确定性减少了。我们可以用韦恩图来可视化这些熵之间的关系from matplotlib_venn import venn2 plt.figure(figsize(6,6)) v venn2(subsets(hx, hy, hxy-hx-hy), set_labels(H(X), H(Y))) plt.title(熵的韦恩图表示) plt.show()3. 互信息量化变量间的依赖关系互信息衡量的是知道一个变量后另一个变量不确定性的减少量。继续上面的天气与穿衣例子# 互信息 I(X;Y) H(Y) - H(Y|X) mutual_info hy - hy_given_x print(f互信息 I(X;Y): {mutual_info:.3f} bits)互信息的值告诉我们天气和穿衣选择之间的关联强度。我们可以修改概率分布观察互信息的变化def plot_mutual_info_heatmap(): p_rain np.linspace(0, 1, 20) p_short_given_rain np.linspace(0, 1, 20) mi_values np.zeros((20, 20)) for i, pr in enumerate(p_rain): for j, psr in enumerate(p_short_given_rain): # 构造联合分布 p_rain_short pr * psr p_rain_long pr * (1-psr) p_sunny_short (1-pr) * 0.8 # 假设晴天穿短袖概率固定 p_sunny_long (1-pr) * 0.2 joint_prob np.array([[p_sunny_short, p_sunny_long], [p_rain_short, p_rain_long]]) # 计算互信息 hxy entropy(joint_prob.flatten()) hx entropy(joint_prob.sum(axis1)) hy entropy(joint_prob.sum(axis0)) mi hy - (hxy - hx) mi_values[i,j] mi plt.figure(figsize(8,6)) plt.imshow(mi_values, extent[0,1,0,1], originlower, cmapviridis) plt.colorbar(label互信息 (bits)) plt.xlabel(P(短袖|雨天)) plt.ylabel(P(雨天)) plt.title(不同概率分布下的互信息变化) plt.show() plot_mutual_info_heatmap()这个热图展示了不同概率分布下互信息的变化帮助我们直观理解变量间的依赖关系。4. 传递熵量化信息流向传递熵是互信息在时间序列上的扩展用于量化一个时间序列对另一个时间序列的信息流向。考虑两个股票价格序列A和B我们想知道A的价格变化是否会影响B未来的价格变化。def transfer_entropy(x, y, delay1, bins10): 计算从x到y的传递熵 # 离散化时间序列 x_discrete np.digitize(x, np.linspace(min(x), max(x), bins)) y_discrete np.digitize(y, np.linspace(min(y), max(y), bins)) # 准备历史数据 y_past y_discrete[:-delay] y_future y_discrete[delay:] x_past x_discrete[:-delay] # 计算各种概率 hist_y_past, _ np.histogram(y_past, binsbins, range(1,bins1)) p_y_past hist_y_past / hist_y_past.sum() hist_joint, _, _ np.histogram2d(y_future, y_past, binsbins, range[[1,bins1],[1,bins1]]) p_joint hist_joint / hist_joint.sum() hist_joint_cond, _, _, _ np.histogramdd(np.column_stack([y_future, y_past, x_past]), binsbins, range[[1,bins1],[1,bins1],[1,bins1]]) p_joint_cond hist_joint_cond / hist_joint_cond.sum() # 计算各项熵 h_y_future_given_y_past entropy(p_joint.flatten()) - entropy(p_y_past) h_y_future_given_y_past_x_past entropy(p_joint_cond.flatten()) - entropy(p_joint_cond.sum(axis(0,1))) # 传递熵 te h_y_future_given_y_past - h_y_future_given_y_past_x_past return max(te, 0) # 传递熵不会为负 # 生成示例数据 np.random.seed(42) n 500 x np.cumsum(np.random.randn(n)) y np.roll(x, 1) * 0.5 np.random.randn(n) * 0.2 # y受x的滞后影响 # 计算传递熵 te_xy transfer_entropy(x, y) te_yx transfer_entropy(y, x) print(fX→Y的传递熵: {te_xy:.3f} bits) print(fY→X的传递熵: {te_yx:.3f} bits)这个例子展示了如何量化两个时间序列间的信息流向。我们可以清楚地看到x对y的影响大于y对x的影响这与我们生成数据的方式一致。5. 实际应用分析脑电信号间的信息流动传递熵在神经科学中有广泛应用可以用来分析不同脑区间的信息流动。下面是一个简化的脑电信号分析示例def simulate_eeg(n_samples1000, n_channels3): 模拟多通道脑电信号 signals np.zeros((n_samples, n_channels)) for i in range(1, n_samples): # 通道0是自发活动 signals[i, 0] 0.8 * signals[i-1, 0] np.random.randn() * 0.5 # 通道1受通道0影响 signals[i, 1] 0.6 * signals[i-1, 1] 0.3 * signals[i-1, 0] np.random.randn() * 0.4 # 通道2受通道1影响 signals[i, 2] 0.7 * signals[i-1, 2] 0.2 * signals[i-1, 1] np.random.randn() * 0.3 return signals eeg_data simulate_eeg() # 计算所有通道间的传递熵 n_channels eeg_data.shape[1] te_matrix np.zeros((n_channels, n_channels)) for i in range(n_channels): for j in range(n_channels): if i ! j: te_matrix[i,j] transfer_entropy(eeg_data[:,i], eeg_data[:,j], bins15) # 可视化传递熵矩阵 plt.figure(figsize(8,6)) plt.imshow(te_matrix, cmaphot) plt.colorbar(label传递熵 (bits)) plt.xticks(range(n_channels), [f通道{i} for i in range(n_channels)]) plt.yticks(range(n_channels), [f通道{i} for i in range(n_channels)]) plt.title(脑电通道间的传递熵矩阵) plt.show()这个模拟展示了传递熵如何揭示脑区间的信息流动方向。在实际应用中这种分析可以帮助研究者理解不同脑区间的功能连接。