import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize

def entropy(p, base=2):
    """计算熵"""
    p = np.array(p, dtype=float)
    p = p[p > 0]
    return -np.sum(p * np.log(p) / np.log(base))

# 验证：对于3取值的随机变量，均匀分布熵最大
print("验证最大熵原理（3个取值）")
print("-" * 50)

n = 3
max_entropy = np.log2(n)
print(f"理论最大熵: {max_entropy:.4f} bits")
print(f"均匀分布 [1/3, 1/3, 1/3]: {entropy([1/3, 1/3, 1/3]):.4f} bits")

# 生成大量随机分布并验证
np.random.seed(42)
max_found = 0
max_dist = None

for _ in range(10000):
    # 从Dirichlet分布采样
    dist = np.random.dirichlet(np.ones(n))
    h = entropy(dist)
    if h > max_found:
        max_found = h
        max_dist = dist

print(f"\n随机搜索最大熵: {max_found:.4f} bits")
print(f"对应分布: {max_dist.round(4)}")
print(f"与均匀分布的差距: {np.abs(max_dist - 1/n).max():.6f}")

# 可视化：二元分布的熵随参数变化
plt.figure(figsize=(10, 6))
p_values = np.linspace(0.001, 0.999, 1000)
h_values = [entropy([p, 1-p]) for p in p_values]

plt.plot(p_values, h_values, 'b-', linewidth=2.5)
plt.axvline(x=0.5, color='r', linestyle='--', alpha=0.7, label='最大值点 p=0.5')
plt.axhline(y=1, color='g', linestyle='--', alpha=0.7, label='最大熵=1 bit')
plt.xlabel('P(X=0)', fontsize=12)
plt.ylabel('熵 H(X) (bits)', fontsize=12)
plt.title('二元随机变量的熵函数', fontsize=14)
plt.grid(True, alpha=0.3)
plt.legend(fontsize=11)
plt.tight_layout()
plt.savefig('max_entropy_binary.png', dpi=150)
plt.show()

print("\n结论：均匀分布确实达到最大熵")

import numpy as np

def conditional_entropy_chain_demo():
    """
    演示条件熵随条件增多而递减
    """
    # 构造一个简单模型
    # X1, X2 是独立的公平硬币
    # Y = X1 XOR X2 （异或）
    
    # 联合分布
    joint = np.zeros((2, 2, 2))
    for x1 in range(2):
        for x2 in range(2):
            y = x1 ^ x2
            joint[x1, x2, y] = 0.25
    
    # 计算各种熵
    # H(Y)
    marginal_y = joint.sum(axis=(0, 1))
    h_y = entropy(marginal_y)
    
    # H(Y|X1)
    # P(Y|X1=0): 当X1=0时，Y=X2，所以P(Y=0|X1=0)=0.5, P(Y=1|X1=0)=0.5
    # P(Y|X1=1): 当X1=1时，Y=1-X2，所以P(Y=0|X1=1)=0.5, P(Y=1|X1=1)=0.5
    h_y_given_x1 = 1.0  # 仍然是均匀分布
    
    # H(Y|X1,X2): 已知X1和X2，Y完全确定
    h_y_given_x1x2 = 0.0
    
    print("条件熵递减演示")
    print("-" * 50)
    print(f"H(Y) = {h_y:.4f} bits")
    print(f"H(Y|X1) = {h_y_given_x1:.4f} bits")
    print(f"H(Y|X1,X2) = {h_y_given_x1x2:.4f} bits")
    print(f"\n递减关系: H(Y) >= H(Y|X1) >= H(Y|X1,X2)? ", end="")
    print(f"{h_y >= h_y_given_x1 >= h_y_given_x1x2}")
    
    # 另一个例子：部分信息
    print("\n\n部分信息示例")
    print("-" * 50)
    
    # X是骰子结果，Y是X的奇偶性
    # P(X=i) = 1/6, Y = X mod 2
    p_x = np.ones(6) / 6
    # P(Y=0|X) = 1 if X in {2,4,6}, else 0
    # P(Y=1|X) = 1 if X in {1,3,5}, else 0
    
    # H(Y)
    p_y = [0.5, 0.5]  # 一半奇数一半偶数
    h_y = entropy(p_y)
    
    # H(Y|X): 已知X，Y完全确定
    h_y_given_x = 0
    
    # H(X|Y): 已知奇偶，还有3种可能
    # P(X=i|Y=0) = 1/3 for i in {2,4,6}
    h_x_given_y_even = entropy([1/3, 1/3, 1/3])
    h_x_given_y = h_x_given_y_even  # 对称
    
    print(f"H(Y) = {h_y:.4f} bits")
    print(f"H(Y|X) = {h_y_given_x:.4f} bits (已知点数，奇偶确定)")
    print(f"H(X|Y) = {h_x_given_y:.4f} bits (已知奇偶，还剩3种可能)")
    print(f"H(X) = {entropy(p_x):.4f} bits")
    print(f"\n互信息 I(X;Y) = H(X) - H(X|Y) = {entropy(p_x) - h_x_given_y:.4f} bits")

conditional_entropy_chain_demo()

import numpy as np
from itertools import product

def verify_chain_rule():
    """验证熵和互信息的链式法则"""
    
    # 构造一个三变量联合分布
    # X, Y 是独立的公平硬币
    # Z = X XOR Y
    
    joint = np.zeros((2, 2, 2))
    for x in range(2):
        for y in range(2):
            z = x ^ y
            joint[x, y, z] = 0.25
    
    # 计算边缘分布
    marginal_xy = joint.sum(axis=2)
    marginal_x = marginal_xy.sum(axis=1)
    marginal_y = marginal_xy.sum(axis=0)
    
    # 计算各种熵
    def calc_entropy(dist):
        flat = dist.flatten()
        flat = flat[flat > 0]
        return -np.sum(flat * np.log2(flat))
    
    h_x = entropy(marginal_x)
    h_y = entropy(marginal_y)
    h_xy = calc_entropy(marginal_xy)
    
    # H(Z)
    marginal_z = joint.sum(axis=(0, 1))
    h_z = entropy(marginal_z)
    
    # H(X,Y,Z)
    h_xyz = calc_entropy(joint)
    
    print("熵的链式法则验证")
    print("-" * 50)
    print(f"H(X) = {h_x:.4f}")
    print(f"H(Y) = {h_y:.4f}")
    print(f"H(Z) = {h_z:.4f}")
    print(f"H(X,Y) = {h_xy:.4f}")
    print(f"H(X,Y,Z) = {h_xyz:.4f}")
    
    # 验证: H(X,Y,Z) = H(X) + H(Y|X) + H(Z|X,Y)
    # H(Y|X) = H(X,Y) - H(X)
    h_y_given_x = h_xy - h_x
    # H(Z|X,Y) = 0 (因为Z由X,Y确定)
    h_z_given_xy = 0
    
    lhs = h_xyz
    rhs = h_x + h_y_given_x + h_z_given_xy
    
    print(f"\n验证: H(X,Y,Z) = H(X) + H(Y|X) + H(Z|X,Y)")
    print(f"左边 = {lhs:.4f}")
    print(f"右边 = {h_x:.4f} + {h_y_given_x:.4f} + {h_z_given_xy:.4f} = {rhs:.4f}")
    print(f"相等? {np.isclose(lhs, rhs)}")
    
    # 互信息的链式法则
    # I(X,Y;Z) = I(X;Z) + I(Y;Z|X)
    
    def calc_mutual_info(joint_2d):
        marginal_x = joint_2d.sum(axis=1)
        marginal_y = joint_2d.sum(axis=0)
        mi = 0
        for i in range(joint_2d.shape[0]):
            for j in range(joint_2d.shape[1]):
                pxy = joint_2d[i, j]
                px = marginal_x[i]
                py = marginal_y[j]
                if pxy > 0 and px > 0 and py > 0:
                    mi += pxy * np.log2(pxy / (px * py))
        return mi
    
    # I(X;Z)
    joint_xz = joint.sum(axis=1)
    i_xz = calc_mutual_info(joint_xz)
    
    # I(Y;Z)
    joint_yz = joint.sum(axis=0)
    i_yz = calc_mutual_info(joint_yz)
    
    # I(X,Y;Z)
    # 需要将(X,Y)视为单个变量
    joint_xyz_flat = joint.reshape(4, 2)  # 4种(X,Y)组合，2种Z
    i_xyz_z = calc_mutual_info(joint_xyz_flat)
    
    print(f"\n互信息计算")
    print("-" * 50)
    print(f"I(X;Z) = {i_xz:.4f}")
    print(f"I(Y;Z) = {i_yz:.4f}")
    print(f"I(X,Y;Z) = {i_xyz_z:.4f}")
    
    # I(Y;Z|X) = I(X,Y;Z) - I(X;Z)
    i_yz_given_x = i_xyz_z - i_xz
    print(f"I(Y;Z|X) = {i_yz_given_x:.4f}")
    
    print(f"\n验证: I(X,Y;Z) = I(X;Z) + I(Y;Z|X)")
    print(f"{i_xyz_z:.4f} = {i_xz:.4f} + {i_yz_given_x:.4f}? {np.isclose(i_xyz_z, i_xz + i_yz_given_x)}")

verify_chain_rule()

import numpy as np

def data_processing_inequality_demo():
    """
    演示数据处理不等式
    
    场景：X是原始信号，Y是通过噪声信道接收的信号，Z是对Y的量化
    """
    
    # X: 4个可能的信号值
    # Y: 接收到的信号（有噪声）
    # Z: 对Y的粗量化
    
    # P(X)
    p_x = np.array([0.25, 0.25, 0.25, 0.25])
    
    # P(Y|X): 噪声信道
    # Y也有4个值，但有噪声
    p_y_given_x = np.array([
        [0.7, 0.1, 0.1, 0.1],  # X=0时，Y大多是0
        [0.1, 0.7, 0.1, 0.1],  # X=1时，Y大多是1
        [0.1, 0.1, 0.7, 0.1],  # ...
        [0.1, 0.1, 0.1, 0.7]
    ])
    
    # P(Z|Y): 量化器（将4个值量化为2个）
    # Z=0 if Y in {0,1}, Z=1 if Y in {2,3}
    p_z_given_y = np.array([
        [1, 0],  # Y=0 -> Z=0
        [1, 0],  # Y=1 -> Z=0
        [0, 1],  # Y=2 -> Z=1
        [0, 1]   # Y=3 -> Z=1
    ])
    
    # 计算联合分布 P(X,Y)
    joint_xy = np.zeros((4, 4))
    for i in range(4):
        for j in range(4):
            joint_xy[i, j] = p_x[i] * p_y_given_x[i, j]
    
    # 计算联合分布 P(Y,Z)
    marginal_y = joint_xy.sum(axis=0)
    joint_yz = np.zeros((4, 2))
    for i in range(4):
        for j in range(2):
            joint_yz[i, j] = marginal_y[i] * p_z_given_y[i, j]
    
    # 计算联合分布 P(X,Z)
    # P(X,Z) = sum_Y P(X,Y)P(Z|Y)
    joint_xz = np.zeros((4, 2))
    for i in range(4):
        for k in range(2):
            for j in range(4):
                joint_xz[i, k] += joint_xy[i, j] * p_z_given_y[j, k]
    
    # 计算互信息
    def mutual_info(joint):
        mx = joint.sum(axis=1)
        my = joint.sum(axis=0)
        mi = 0
        for i in range(joint.shape[0]):
            for j in range(joint.shape[1]):
                pxy = joint[i, j]
                px = mx[i]
                py = my[j]
                if pxy > 0 and px > 0 and py > 0:
                    mi += pxy * np.log2(pxy / (px * py))
        return mi
    
    i_xy = mutual_info(joint_xy)
    i_xz = mutual_info(joint_xz)
    
    print("数据处理不等式演示")
    print("-" * 50)
    print("场景: X(原始) -> Y(噪声) -> Z(量化)")
    print(f"\nI(X;Y) = {i_xy:.4f} bits")
    print(f"I(X;Z) = {i_xz:.4f} bits")
    print(f"\nI(X;Y) >= I(X;Z)? {i_xy >= i_xz}")
    print(f"信息损失: {i_xy - i_xz:.4f} bits")
    
    # 改变噪声水平
    print("\n\n不同噪声水平下的信息损失")
    print("-" * 50)
    
    noise_levels = [0.01, 0.1, 0.2, 0.3, 0.4]
    for noise in noise_levels:
        # 调整噪声信道
        p_ygx = np.array([
            [1-3*noise, noise, noise, noise],
            [noise, 1-3*noise, noise, noise],
            [noise, noise, 1-3*noise, noise],
            [noise, noise, noise, 1-3*noise]
        ])
        
        joint_xy = np.zeros((4, 4))
        for i in range(4):
            for j in range(4):
                joint_xy[i, j] = p_x[i] * p_ygx[i, j]
        
        marginal_y = joint_xy.sum(axis=0)
        joint_xz = np.zeros((4, 2))
        for i in range(4):
            for k in range(2):
                for j in range(4):
                    joint_xz[i, k] += joint_xy[i, j] * p_z_given_y[j, k]
        
        i_xy = mutual_info(joint_xy)
        i_xz = mutual_info(joint_xz)
        
        print(f"噪声={noise:.2f}: I(X;Y)={i_xy:.4f}, I(X;Z)={i_xz:.4f}, 损失={i_xy-i_xz:.4f}")

data_processing_inequality_demo()

import numpy as np

def fano_inequality_demo():
    """
    验证Fano不等式
    """
    
    # 场景：X是发送的符号，Y是接收的信号
    # 基于Y估计X，计算错误概率和条件熵
    
    # 4个符号的通信系统
    n_symbols = 4
    
    # 转移概率矩阵 P(Y|X)
    # 对角线是正确接收的概率
    correct_prob = 0.8
    error_prob = (1 - correct_prob) / (n_symbols - 1)
    
    p_y_given_x = np.full((n_symbols, n_symbols), error_prob)
    for i in range(n_symbols):
        p_y_given_x[i, i] = correct_prob
    
    # 均匀先验
    p_x = np.ones(n_symbols) / n_symbols
    
    # 联合分布
    joint = np.zeros((n_symbols, n_symbols))
    for i in range(n_symbols):
        for j in range(n_symbols):
            joint[i, j] = p_x[i] * p_y_given_x[i, j]
    
    # MAP估计（最大后验）
    # 由于对称性，MAP就是选择最可能的Y对应的X
    estimates = np.argmax(joint, axis=0)  # 对每个Y，选择最可能的X
    
    # 计算错误概率
    p_error = 0
    for j in range(n_symbols):
        p_y = joint[:, j].sum()
        for i in range(n_symbols):
            if i != estimates[j]:
                p_error += joint[i, j]
    
    # 计算条件熵 H(X|Y)
    def conditional_entropy_from_joint(joint):
        n_x, n_y = joint.shape
        marginal_y = joint.sum(axis=0)
        
        h = 0
        for j in range(n_y):
            if marginal_y[j] > 0:
                cond_dist = joint[:, j] / marginal_y[j]
                cond_dist = cond_dist[cond_dist > 0]
                h_y_given_j = -np.sum(cond_dist * np.log2(cond_dist))
                h += marginal_y[j] * h_y_given_j
        
        return h
    
    h_x_given_y = conditional_entropy_from_joint(joint)
    
    # Fano不等式的右边
    h_pe = -p_error * np.log2(p_error + 1e-10) - (1-p_error) * np.log2(1-p_error + 1e-10)
    fanos_rhs = h_pe + p_error * np.log2(n_symbols - 1)
    
    print("Fano不等式验证")
    print("=" * 60)
    print(f"符号数: {n_symbols}")
    print(f"正确接收概率: {correct_prob}")
    print(f"\n错误概率 P_e = {p_error:.4f}")
    print(f"条件熵 H(X|Y) = {h_x_given_y:.4f} bits")
    print(f"\nFano不等式右边:")
    print(f"  H(P_e) = {h_pe:.4f}")
    print(f"  P_e * log2(|X|-1) = {p_error:.4f} * {np.log2(n_symbols-1):.4f} = {p_error * np.log2(n_symbols-1):.4f}")
    print(f"  总和 = {fanos_rhs:.4f}")
    print(f"\nH(X|Y) <= Fano上界? {h_x_given_y <= fanos_rhs + 1e-6}")
    print(f"差距: {fanos_rhs - h_x_given_y:.4f} bits")
    
    # 测试不同正确概率
    print("\n\n不同信道质量下的Fano不等式")
    print("=" * 60)
    
    for cp in [0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 0.99]:
        ep = (1 - cp) / (n_symbols - 1)
        p_ygx = np.full((n_symbols, n_symbols), ep)
        for i in range(n_symbols):
            p_ygx[i, i] = cp
        
        joint = np.zeros((n_symbols, n_symbols))
        for i in range(n_symbols):
            for j in range(n_symbols):
                joint[i, j] = p_x[i] * p_ygx[i, j]
        
        pe = 0
        for j in range(n_symbols):
            for i in range(n_symbols):
                if i != np.argmax(joint[:, j]):
                    pe += joint[i, j]
        
        hxgy = conditional_entropy_from_joint(joint)
        hpe = -pe * np.log2(pe + 1e-10) - (1-pe) * np.log2(1-pe + 1e-10)
        fanos = hpe + pe * np.log2(n_symbols - 1)
        
        print(f"P(correct)={cp:.2f}: P_e={pe:.4f}, H(X|Y)={hxgy:.4f}, Fano={fanos:.4f}, 满足={hxgy <= fanos + 1e-6}")

fano_inequality_demo()

import numpy as np
import matplotlib.pyplot as plt

def verify_entropy_concavity():
    """验证熵函数的凹性"""
    
    # 两个极端分布
    p1 = np.array([0.9, 0.05, 0.05])
    p2 = np.array([0.05, 0.05, 0.9])
    
    h1 = entropy(p1)
    h2 = entropy(p2)
    
    print("熵的凹性验证")
    print("=" * 60)
    print(f"分布1: {p1}, H = {h1:.4f}")
    print(f"分布2: {p2}, H = {h2:.4f}")
    
    # 测试不同的混合比例
    lambdas = np.linspace(0, 1, 11)
    h_mixed = []
    h_weighted = []
    
    for lam in lambdas:
        p_mix = lam * p1 + (1 - lam) * p2
        h_m = entropy(p_mix)
        h_w = lam * h1 + (1 - lam) * h2
        
        h_mixed.append(h_m)
        h_weighted.append(h_w)
        
        print(f"λ={lam:.1f}: H(混合)={h_m:.4f}, λH₁+(1-λ)H₂={h_w:.4f}, 差={h_m-h_w:.4f}")
    
    # 可视化
    plt.figure(figsize=(10, 6))
    plt.plot(lambdas, h_mixed, 'bo-', linewidth=2, markersize=8, label='H(λp₁ + (1-λ)p₂)')
    plt.plot(lambdas, h_weighted, 'r--', linewidth=2, label='λH(p₁) + (1-λ)H(p₂)')
    plt.xlabel('λ', fontsize=12)
    plt.ylabel('熵 (bits)', fontsize=12)
    plt.title('熵函数的凹性', fontsize=14)
    plt.legend(fontsize=11)
    plt.grid(True, alpha=0.3)
    plt.fill_between(lambdas, h_weighted, h_mixed, alpha=0.2, color='green', 
                      label='凹性间隙')
    plt.tight_layout()
    plt.savefig('entropy_concavity.png', dpi=150)
    plt.show()
    
    print(f"\n所有点都满足 H(混合) ≥ 加权平均? {all(m >= w - 1e-10 for m, w in zip(h_mixed, h_weighted))}")

verify_entropy_concavity()

import numpy as np
from sklearn.linear_model import LogisticRegression
from sklearn.datasets import make_classification

# 生成数据
X, y = make_classification(n_samples=1000, n_features=2, n_redundant=0, 
                           random_state=42)

# 训练逻辑回归（使用交叉熵损失）
model = LogisticRegression()
model.fit(X, y)

# 计算训练集交叉熵
y_prob = model.predict_proba(X)
cross_entropy = -np.mean(y * np.log(y_prob[:, 1] + 1e-10) + 
                         (1-y) * np.log(y_prob[:, 0] + 1e-10))

print(f"交叉熵损失: {cross_entropy:.4f}")
print(f"准确率: {model.score(X, y):.4f}")

from sklearn.feature_selection import mutual_info_classif
from sklearn.datasets import load_iris

iris = load_iris()
X, y = iris.data, iris.target

# 计算每个特征与标签的互信息
mi_scores = mutual_info_classif(X, y, random_state=42)

print("Iris数据集特征互信息:")
for name, score in zip(iris.feature_names, mi_scores):
    print(f"  {name:20s}: {score:.4f}")

# 基于互信息的特征选择
selected_features = np.argsort(mi_scores)[-2:]  # 选择前2个
print(f"\n选择的特征索引: {selected_features}")

import numpy as np

class InformationInequalityVerifier:
    """信息论不等式验证器"""
    
    @staticmethod
    def verify_non_negativity(p):
        """验证熵的非负性"""
        h = entropy(p)
        return h >= 0, h
    
    @staticmethod
    def verify_max_entropy(p):
        """验证最大熵原理"""
        h = entropy(p)
        h_max = np.log2(len(p))
        return h <= h_max + 1e-10, h, h_max
    
    @staticmethod
    def verify_gibbs(p, q):
        """验证吉布斯不等式 D_KL(p||q) >= 0"""
        kl = kl_divergence(p, q)
        return kl >= -1e-10, kl
    
    @staticmethod
    def verify_chain_rule(joint):
        """验证熵的链式法则"""
        joint = np.array(joint)
        
        # H(X,Y)
        h_xy = calc_joint_entropy(joint)
        
        # H(X)
        marginal_x = joint.sum(axis=1)
        h_x = entropy(marginal_x)
        
        # H(Y|X)
        h_y_given_x = calc_conditional_entropy(joint, axis=0)
        
        lhs = h_xy
        rhs = h_x + h_y_given_x
        
        return np.isclose(lhs, rhs), lhs, rhs
    
    @staticmethod
    def verify_data_processing(joint_xyz):
        """验证数据处理不等式"""
        # 实现略
        pass


def calc_joint_entropy(joint):
    flat = np.array(joint).flatten()
    flat = flat[flat > 0]
    return -np.sum(flat * np.log2(flat))

def calc_conditional_entropy(joint, axis=0):
    joint = np.array(joint)
    if axis == 0:
        marginal = joint.sum(axis=1)
    else:
        marginal = joint.sum(axis=0)
    
    h = 0
    for i in range(joint.shape[axis]):
        if axis == 0:
            cond = joint[i, :]
        else:
            cond = joint[:, i]
        
        px = marginal[i]
        if px > 0:
            cond_prob = cond / px
            cond_prob = cond_prob[cond_prob > 0]
            h += px * (-np.sum(cond_prob * np.log2(cond_prob)))
    
    return h

def kl_divergence(p, q):
    p = np.array(p)
    q = np.array(q)
    mask = p > 0
    return np.sum(p[mask] * np.log2(p[mask] / q[mask]))


# 测试
verifier = InformationInequalityVerifier()

print("非负性验证:", verifier.verify_non_negativity([0.5, 0.3, 0.2]))
print("最大熵验证:", verifier.verify_max_entropy([0.5, 0.3, 0.2]))
print("吉布斯不等式:", verifier.verify_gibbs([0.5, 0.5], [0.3, 0.7]))

import numpy as np
import matplotlib.pyplot as plt

class MarkovChainSimulator:
    """马尔可夫链模拟器，用于演示数据处理不等式"""
    
    def __init__(self, transition_matrix):
        self.P = np.array(transition_matrix)
        self.n_states = self.P.shape[0]
    
    def simulate(self, initial_state, n_steps):
        """模拟马尔可夫链"""
        states = [initial_state]
        current = initial_state
        
        for _ in range(n_steps - 1):
            current = np.random.choice(self.n_states, p=self.P[current])
            states.append(current)
        
        return states
    
    def compute_stationary_distribution(self):
        """计算平稳分布"""
        eigenvalues, eigenvectors = np.linalg.eig(self.P.T)
        # 找到对应于特征值1的特征向量
        idx = np.argmin(np.abs(eigenvalues - 1))
        pi = np.real(eigenvectors[:, idx])
        pi = pi / pi.sum()
        return pi
    
    def compute_mutual_information_over_time(self, initial_dist):
        """计算I(X_0; X_t)随时间的衰减"""
        pi = np.array(initial_dist)
        mi_values = []
        
        # P^t
        P_power = np.eye(self.n_states)
        
        for t in range(20):
            # P(X_t | X_0) = P^t
            conditional = P_power
            
            # 联合分布 P(X_0, X_t)
            joint = np.zeros((self.n_states, self.n_states))
            for i in range(self.n_states):
                for j in range(self.n_states):
                    joint[i, j] = pi[i] * conditional[i, j]
            
            # 计算互信息
            mi = self._mutual_info(joint)
            mi_values.append(mi)
            
            P_power = P_power @ self.P
        
        return mi_values
    
    def _mutual_info(self, joint):
        mx = joint.sum(axis=1)
        my = joint.sum(axis=0)
        mi = 0
        for i in range(joint.shape[0]):
            for j in range(joint.shape[1]):
                pxy = joint[i, j]
                px = mx[i]
                py = my[j]
                if pxy > 0 and px > 0 and py > 0:
                    mi += pxy * np.log2(pxy / (px * py))
        return mi


# 示例：三状态马尔可夫链
P = np.array([
    [0.7, 0.2, 0.1],
    [0.1, 0.8, 0.1],
    [0.1, 0.1, 0.8]
])

simulator = MarkovChainSimulator(P)
initial = [0.5, 0.3, 0.2]

mi_over_time = simulator.compute_mutual_information_over_time(initial)

plt.figure(figsize=(10, 6))
plt.plot(mi_over_time, 'bo-', linewidth=2, markersize=8)
plt.xlabel('时间步 t', fontsize=12)
plt.ylabel('I(X₀; Xₜ) (bits)', fontsize=12)
plt.title('互信息随时间的衰减（数据处理不等式）', fontsize=14)
plt.grid(True, alpha=0.3)
plt.axhline(y=0, color='r', linestyle='--', alpha=0.5)
plt.tight_layout()
plt.savefig('mi_decay.png', dpi=150)
plt.show()

print("平稳分布:", simulator.compute_stationary_distribution())
print("互信息最终趋于0，体现了数据处理不等式")

import numpy as np
import matplotlib.pyplot as plt

def visualize_fano_bound(n_symbols=4):
    """可视化Fano不等式"""
    
    p_errors = np.linspace(0.01, 0.99, 100)
    
    # Fano上界
    fanos_bounds = []
    for pe in p_errors:
        h_pe = -pe * np.log2(pe) - (1-pe) * np.log2(1-pe)
        bound = h_pe + pe * np.log2(n_symbols - 1)
        fanos_bounds.append(bound)
    
    plt.figure(figsize=(10, 6))
    plt.plot(p_errors, fanos_bounds, 'b-', linewidth=2.5, label='Fano上界')
    plt.xlabel('错误概率 P_e', fontsize=12)
    plt.ylabel('H(X|Ŷ) 上界 (bits)', fontsize=12)
    plt.title(f'Fano不等式 ({n_symbols}个符号)', fontsize=14)
    plt.grid(True, alpha=0.3)
    plt.legend(fontsize=11)
    plt.tight_layout()
    plt.savefig('fano_bound.png', dpi=150)
    plt.show()

visualize_fano_bound(4)