02 - 随机过程与马尔可夫链¶

import numpy as np
import matplotlib.pyplot as plt

# 定义天气马尔可夫链
states = ['Sunny', 'Rainy']
P = np.array([  # np.array创建NumPy数组
    [0.9, 0.1],  # 从Sunny出发
    [0.5, 0.5]   # 从Rainy出发
])

def simulate_markov_chain(P, initial_state, n_steps):
    """
    模拟马尔可夫链

    参数:
        P: 转移矩阵
        initial_state: 初始状态索引
        n_steps: 模拟步数
    """
    n_states = P.shape[0]
    states_history = [initial_state]

    current_state = initial_state
    for _ in range(n_steps):
        # 根据当前状态的概率分布采样下一个状态
        next_state = np.random.choice(n_states, p=P[current_state])
        states_history.append(next_state)
        current_state = next_state

    return states_history

# 模拟
np.random.seed(42)
n_steps = 100
initial = 0  # 从晴天开始

history = simulate_markov_chain(P, initial, n_steps)

# 可视化
plt.figure(figsize=(12, 4))

plt.subplot(1, 2, 1)
plt.plot(history, 'b-', linewidth=1)
plt.yticks([0, 1], states)
plt.xlabel('Time Step')
plt.ylabel('State')
plt.title('Markov Chain Simulation')
plt.grid(True, alpha=0.3)

# 统计各状态出现频率
plt.subplot(1, 2, 2)
unique, counts = np.unique(history, return_counts=True)
plt.bar([states[i] for i in unique], counts / len(history))
plt.ylabel('Frequency')
plt.title('State Distribution')
plt.grid(True, alpha=0.3, axis='y')

plt.tight_layout()
plt.savefig('markov_chain_weather.png', dpi=150)
plt.show()

print(f"初始状态: {states[initial]}")
print(f"晴天频率: {counts[0]/len(history):.2%}")
print(f"雨天频率: {counts[1]/len(history):.2%}")

# 方法1：模拟法（大数定律）
def estimate_stationary_distribution(P, n_steps=10000):
    """通过长时间模拟估计稳态分布"""
    n_states = P.shape[0]
    state_counts = np.zeros(n_states)

    current_state = 0
    for _ in range(n_steps):
        state_counts[current_state] += 1
        current_state = np.random.choice(n_states, p=P[current_state])

    return state_counts / n_steps

# 方法2：特征值法（精确解）
def compute_stationary_distribution(P):
    """通过求解特征向量计算稳态分布"""
    # 求解 πP = π，即 P^T π^T = π^T
    eigenvalues, eigenvectors = np.linalg.eig(P.T)  # np.linalg线性代数运算

    # 找到特征值为1的特征向量
    idx = np.argmin(np.abs(eigenvalues - 1))
    stationary = eigenvectors[:, idx].real

    # 归一化
    stationary = stationary / np.sum(stationary)
    return stationary

# 计算
stationary_sim = estimate_stationary_distribution(P)
stationary_exact = compute_stationary_distribution(P)

print("稳态分布:")
print(f"  模拟估计: {stationary_sim}")
print(f"  精确计算: {stationary_exact}")
print(f"  验证: πP = {stationary_exact @ P}")

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

def forward_diffusion_1d(x_0, n_steps, beta_schedule='linear'):
    """
    一维扩散过程模拟

    参数:
        x_0: 初始值
        n_steps: 扩散步数
        beta_schedule: beta调度策略
    """
    if beta_schedule == 'linear':
        betas = np.linspace(0.01, 0.99, n_steps)
    else:
        betas = np.ones(n_steps) * 0.1

    # 计算alpha
    alphas = 1 - betas
    alphas_cumprod = np.cumprod(alphas)

    # 存储每一步的状态
    trajectory = [x_0]
    x = x_0

    for t in range(n_steps):
        # x_t = sqrt(1-beta_t) * x_{t-1} + sqrt(beta_t) * epsilon
        epsilon = np.random.normal(0, 1)
        x = np.sqrt(1 - betas[t]) * x + np.sqrt(betas[t]) * epsilon
        trajectory.append(x)

    return np.array(trajectory), alphas_cumprod

# 模拟多个轨迹
np.random.seed(42)
n_steps = 100
n_trajectories = 10

plt.figure(figsize=(12, 5))

plt.subplot(1, 2, 1)
for i in range(n_trajectories):
    x_0 = np.random.normal(0, 1)
    traj, _ = forward_diffusion_1d(x_0, n_steps)
    plt.plot(traj, alpha=0.6, linewidth=1)

plt.xlabel('Time Step')
plt.ylabel('x value')
plt.title('Forward Diffusion Process (10 trajectories)')
plt.grid(True, alpha=0.3)
plt.axhline(y=0, color='r', linestyle='--', alpha=0.5, label='Mean')

# 显示分布随时间的变化
plt.subplot(1, 2, 2)
n_samples = 1000
samples_t = {0: [], 25: [], 50: [], 75: [], 99: []}

for _ in range(n_samples):
    x_0 = np.random.normal(0, 1)
    traj, alphas_cumprod = forward_diffusion_1d(x_0, n_steps)
    for t in samples_t.keys():
        samples_t[t].append(traj[t])

# 绘制不同时间步的分布
for i, (t, samples) in enumerate(samples_t.items()):  # enumerate同时获取索引和元素
    plt.hist(samples, bins=30, alpha=0.5,
             label=f't={t}, α={alphas_cumprod[t]:.3f}')

plt.xlabel('x value')
plt.ylabel('Count')
plt.title('Distribution at Different Time Steps')
plt.legend()
plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('forward_diffusion_1d.png', dpi=150)
plt.show()

print("观察：")
print("1. 随着时间步增加，所有轨迹趋向于0附近")
print("2. 分布逐渐变成标准高斯分布")
print("3. 这就是马尔可夫链收敛到稳态分布的过程")

def q_sample_direct(x_0, t, alphas_cumprod):
    """
    直接从 q(x_t | x_0) 采样
    利用马尔可夫链的性质

    参数:
        x_0: 初始值
        t: 目标时间步
        alphas_cumprod: 累积alpha值
    """
    alpha_bar_t = alphas_cumprod[t]

    # 均值: sqrt(alpha_bar_t) * x_0
    mean = np.sqrt(alpha_bar_t) * x_0

    # 标准差: sqrt(1 - alpha_bar_t)
    std = np.sqrt(1 - alpha_bar_t)

    # 采样
    epsilon = np.random.normal(0, 1, x_0.shape)
    x_t = mean + std * epsilon

    return x_t, epsilon

# 验证：直接采样 vs 逐步模拟
np.random.seed(42)
n_steps = 100
betas = np.linspace(0.01, 0.99, n_steps)
alphas = 1 - betas
alphas_cumprod = np.cumprod(alphas)

x_0 = np.array([2.0])
t = 50

# 方法1：逐步模拟
x_step = x_0.copy()
for i in range(t):
    epsilon = np.random.normal(0, 1, x_0.shape)
    x_step = np.sqrt(alphas[i]) * x_step + np.sqrt(betas[i]) * epsilon

# 方法2：直接采样
np.random.seed(42)  # 重置种子，使用相同的随机数
x_direct, _ = q_sample_direct(x_0, t, alphas_cumprod)

print(f"初始值 x_0 = {x_0[0]}")
print(f"逐步模拟 x_{t} = {x_step[0]:.4f}")
print(f"直接采样 x_{t} = {x_direct[0]:.4f}")
print(f"两者相等！这证明了马尔可夫链的性质")