""" ============================================================================= 第 1 讲 · 代码 2：马尔可夫链（Markov Chain / Markov Process） ============================================================================= 对应讲义：lectures/lecture01_intro_to_rl.md 第 3.5 节（马尔可夫假设）、第 7 节（马尔可夫过程） ----------------------------------------------------------------------------- 【这个文件要解决什么问题？】 ----------------------------------------------------------------------------- 马尔可夫链是 RL 数学的第一块积木：只有"状态"和"状态如何随机转移"， 没有动作、没有奖励。我们用火星车的例子，亲手完成下面几件事： 1. 用一个 7×7 的矩阵表示转移概率，并检查它是否合法 2. 按转移概率"掷骰子"，采样出一条条轨迹 3. 计算"t 步之后火星车在各个位置的概率"，观察它收敛到平稳分布 4. 只给数据、不给矩阵，通过"数数"把转移矩阵估计出来（学习模型！） 5. 用数据验证马尔可夫性："上一步在哪"不影响"下一步去哪" ----------------------------------------------------------------------------- 【为什么要学这个？】 ----------------------------------------------------------------------------- 当策略固定时，MDP 就退化成马尔可夫链（加上奖励就是 MRP）。 所以"采样轨迹""状态分布""估计转移模型"这些操作，在后面每一讲都会反复出现： - 采样轨迹 → 第 3 讲蒙特卡洛方法的数据来源 - 状态分布 → 第 5、6 讲函数逼近、策略梯度中"按访问频率加权" - 估计转移模型 → 第 3 讲确定性等价方法、第 13 讲基于模型的探索、第 15 讲 Dyna ----------------------------------------------------------------------------- 【运行方式】 ----------------------------------------------------------------------------- python 02_markov_chain.py 【下标约定】代码中的状态 i <==> 讲义中的 s_{i+1} """ import numpy as np # 【Python 小知识】打印 numpy 数组时的显示设置： # precision=3 → 保留 3 位小数；suppress=True → 不用科学计数法（比如不显示 1e-05）。 # 这只影响"打印出来长什么样"，不影响计算精度。 np.set_printoptions(precision=3, suppress=True) # ============================================================================= # 第 1 部分：构建并检查转移矩阵 # ============================================================================= def build_mars_rover_chain(): """ 构建讲义 7.4 节的火星车马尔可夫链转移矩阵。 【思路】 转移矩阵 P 是一个 N×N 的方阵，P[i, j] = "现在在 i，下一步到 j 的概率"。 行 = 从哪里来，列 = 到哪里去。 火星车链的规则： - 最左格 0： 0.6 留在原地，0.4 向右 - 中间格 1~5：0.4 向左，0.2 原地，0.4 向右 - 最右格 6： 0.6 留在原地，0.4 向左 【返回】 P : numpy 数组，形状 (7, 7)，每行加起来等于 1。 【Python 小知识】 np.zeros((7, 7)) 创建一个 7 行 7 列、全是 0 的矩阵。 注意括号：形状 (7, 7) 本身是一个元组，所以要写两层括号。 P[i, j] = x 把第 i 行第 j 列设为 x（行列都从 0 开始数）。 """ n = 7 P = np.zeros((n, n)) # 最左边的格子 P[0, 0] = 0.6 P[0, 1] = 0.4 # 中间的格子 1, 2, 3, 4, 5 # 【Python 小知识】range(1, n - 1) 生成 1, 2, ..., n-2（不含 n-1），即 1~5。 for i in range(1, n - 1): P[i, i - 1] = 0.4 # 向左 P[i, i] = 0.2 # 原地 P[i, i + 1] = 0.4 # 向右 # 最右边的格子 P[n - 1, n - 1] = 0.6 P[n - 1, n - 2] = 0.4 return P def check_transition_matrix(P, atol=1e-8): """ 检查一个矩阵是否是合法的转移矩阵（随机矩阵）。不合法就报错，合法就什么都不返回。 【思路】 合法的转移矩阵必须满足三个条件（讲义 7.3 节）： (1) 是方阵（N×N）——"从 N 个状态出发，到 N 个状态去" (2) 所有元素 >= 0——概率不能是负数 (3) 每一行之和 = 1——从任何状态出发，下一步总得去某个地方 在写任何用到 P 的算法之前先调用它，能在最早的时候发现"矩阵写错了"这种 bug， 否则你可能会在很久之后看到一个莫名其妙的价值函数，却不知道问题在哪。 【参数】 P : 二维 numpy 数组。 atol : float，允许的浮点误差（absolute tolerance），默认 1e-8。 【边界】 - 为什么需要 atol？因为计算机里的小数是近似存储的。 例如 0.1 + 0.2 在计算机中等于 0.30000000000000004，而不是精确的 0.3。 所以"行和等于 1"必须写成"行和与 1 的差距小于一个很小的数"。 永远不要用 == 比较两个浮点数！ - 空矩阵（0×0）也会被拒绝，因为没有状态的马尔可夫链没有意义。 - 如果 P 里有 NaN（"不是一个数"，通常来自 0/0），也会被拒绝。 【返回】 None（检查通过时）。不通过时抛出 ValueError。 """ # 【Python 小知识】np.asarray(P) 确保 P 是 numpy 数组（如果传入的是嵌套列表，会被转换）。 P = np.asarray(P, dtype=float) # P.ndim 是数组的维数；P.shape 是各维的长度，二维时为 (行数, 列数)。 if P.ndim != 2 or P.shape[0] != P.shape[1]: raise ValueError(f"转移矩阵必须是方阵，但形状是 {P.shape}") if P.shape[0] == 0: raise ValueError("转移矩阵不能为空（至少要有 1 个状态）") # np.isnan(P) 返回一个同形状的布尔矩阵；.any() 判断"是否至少有一个 True"。 if np.isnan(P).any(): raise ValueError("转移矩阵中含有 NaN") if (P < -atol).any(): raise ValueError("转移矩阵中有负数，概率不能为负") # P.sum(axis=1) 对每一行求和，得到一个长度为 N 的向量。 # axis=1 表示"沿着列的方向加"（把同一行的各列加起来）→ 每行一个和； # axis=0 则是把同一列的各行加起来 → 每列一个和。初学者很容易搞反，记住： # axis=1 → 结果长度 = 行数。 row_sums = P.sum(axis=1) # np.allclose(a, b, atol=...) 判断两个数组是否"在误差范围内相等"。 if not np.allclose(row_sums, 1.0, atol=atol): # np.flatnonzero 返回不满足条件的下标，帮助定位是哪一行出了问题。 bad_rows = np.flatnonzero(~np.isclose(row_sums, 1.0, atol=atol)) raise ValueError(f"以下行的和不等于 1：{bad_rows.tolist()}，对应行和为 {row_sums[bad_rows]}") # ============================================================================= # 第 2 部分：采样轨迹 # ============================================================================= def sample_next_state(P, current_state, rng): """ 按照转移概率，从当前状态随机抽取下一个状态。 【思路】 P 的第 current_state 行 P[current_state] 就是"下一步去各个状态的概率"， 是一个概率分布。我们按这个分布抽一个状态。 例：在格子 3，P[3] = [0, 0, 0.4, 0.2, 0.4, 0, 0]， 那么有 40% 概率抽到 2，20% 抽到 3，40% 抽到 4。 【参数】 P : 转移矩阵，形状 (N, N)，应当已通过 check_transition_matrix。 current_state : int，当前状态，必须在 0 到 N-1 之间。 rng : numpy 随机数生成器（np.random.default_rng(...) 创建）。 把 rng 作为参数传进来（而不是在函数内部新建）， 是为了让调用者能控制随机种子，从而复现结果。 【返回】 int，下一个状态。 【边界】 - current_state 越界 → IndexError（由我们主动检查并给出清楚的提示）。 - 本函数为了速度不重复检查 P 的合法性，调用者应事先检查一次。 """ n = P.shape[0] if not (0 <= current_state < n): raise IndexError(f"current_state 必须在 0 到 {n - 1} 之间，但收到了 {current_state}") # 【Python 小知识】rng.choice(n, p=概率向量) 从 0, 1, ..., n-1 中按给定概率抽一个数。 # p 的长度必须等于 n，且元素和为 1（否则 numpy 会报错）。 return int(rng.choice(n, p=P[current_state])) def sample_trajectory(P, start_state, n_steps, rng): """ 从 start_state 出发，连续采样 n_steps 次，得到一条轨迹。 【参数】 P : 转移矩阵，形状 (N, N)。 start_state : int，起点。 n_steps : int，转移的次数，必须 >= 0。 rng : numpy 随机数生成器。 【返回】 list[int]，长度为 n_steps + 1 的状态列表 [s_0, s_1, ..., s_{n_steps}]。 注意长度是 n_steps + 1 而不是 n_steps，因为包含了起点。 （"走 3 步"会经过 4 个位置：起点 + 3 个新位置。这是一个经典的"差一错误"来源。） 【边界】 - n_steps = 0 → 返回 [start_state]，只有起点。 - n_steps < 0 → ValueError。 """ if n_steps < 0: raise ValueError(f"n_steps 必须 >= 0，但收到了 {n_steps}") trajectory = [start_state] state = start_state for _ in range(n_steps): state = sample_next_state(P, state, rng) trajectory.append(state) return trajectory def format_trajectory(trajectory): """ 把状态列表 [3, 4, 5] 转成讲义风格的字符串 "s4 → s5 → s6"，方便阅读。 【Python 小知识】 " → ".join([...]) 用 " → " 把列表里的字符串连接起来。 列表推导式里的 f"s{s + 1}" 把代码编号（从 0 开始）转为讲义编号（从 1 开始）。 """ return " → ".join(f"s{s + 1}" for s in trajectory) # ============================================================================= # 第 3 部分：状态分布的演化与平稳分布 # ============================================================================= def propagate_distribution(P, initial_distribution, n_steps): """ 计算状态分布随时间的演化：d_{t+1} = d_t P。 【思路】 不再追踪"某一条"轨迹，而是追踪"所有可能轨迹"的整体： d_t[s] = "在时刻 t 处于状态 s 的概率"。 由全概率公式：d_{t+1}[s'] = Σ_s d_t[s] * P[s, s']， 写成矩阵乘法就是 d_{t+1} = d_t @ P（行向量乘矩阵）。 【参数】 P : 转移矩阵，形状 (N, N)。 initial_distribution : 长度为 N 的一维数组，初始分布 d_0，元素非负且和为 1。 例：确定从格子 3 出发 → [0, 0, 0, 1, 0, 0, 0]。 n_steps : int，要演化多少步，必须 >= 0。 【返回】 numpy 数组，形状 (n_steps + 1, N)。第 t 行就是 d_t。 【边界】 - initial_distribution 长度不对或不是概率分布 → ValueError。 - n_steps = 0 → 只返回 d_0 一行。 【Python 小知识】 @ 是 Python 的矩阵乘法运算符。对于一维数组 d（长度 N）和二维数组 P（N×N）， d @ P 的结果是长度为 N 的一维数组，第 j 个元素 = Σ_i d[i] * P[i, j]。 注意 @ 和 * 完全不同：* 是"逐元素相乘"，@ 才是线性代数里的矩阵乘法。 """ d = np.asarray(initial_distribution, dtype=float) n = P.shape[0] if d.shape != (n,): raise ValueError(f"initial_distribution 的长度必须是 {n}，但形状是 {d.shape}") if (d < 0).any() or not np.isclose(d.sum(), 1.0): raise ValueError("initial_distribution 必须是概率分布（非负且和为 1）") if n_steps < 0: raise ValueError(f"n_steps 必须 >= 0，但收到了 {n_steps}") history = [d] for _ in range(n_steps): d = d @ P history.append(d) # np.array(列表) 把"由一维数组组成的列表"堆叠成一个二维数组，每个一维数组成为一行。 return np.array(history) def stationary_distribution(P, tol=1e-12, max_iters=100_000): """ 用"幂迭代"求平稳分布 d*（满足 d* = d* P 的分布）。 【思路】 最直接的想法：从任意分布（这里用均匀分布）出发，反复乘 P， 直到分布不再变化。这就是讲义 7.6 节表格里看到的现象： 不管从哪里出发，分布最终都收敛到同一个 d*。 【参数】 tol : float，收敛阈值。相邻两次分布的最大差值小于 tol 就认为收敛。 max_iters : int，最多迭代多少次，防止不收敛时死循环。 【Python 小知识】100_000 就是 100000，下划线只是为了方便人眼读数字。 【返回】 (d, converged) 二元组： - d : 长度为 N 的数组，最后一次迭代得到的分布 - converged : bool，是否在 max_iters 次内收敛 【边界】（重要！） 幂迭代**不是对所有马尔可夫链都收敛**： - 周期链：例如 P = [[0, 1], [1, 0]]（两个状态来回跳）。 如果从 [1, 0] 出发，分布会在 [1,0] 和 [0,1] 之间永远振荡，不会收敛。 （本函数从均匀分布 [0.5, 0.5] 出发，恰好就是它的平稳分布，所以这个例子能"侥幸"收敛； 但对别的周期链，均匀分布不一定是平稳分布，从它出发同样可能永远振荡。） - 可约链：例如有两个互不相通的"区域"，平稳分布不唯一，结果取决于初始分布。 所以返回值里带上 converged 标志，调用者必须检查它。 火星车链是"不可约且非周期"的（任何格子都能互相到达，且有"原地不动"的概率）， 所以保证收敛到唯一的平稳分布。 """ n = P.shape[0] d = np.ones(n) / n # 均匀分布作为起点：np.ones(n) 是全 1 向量，除以 n 后和为 1 for _ in range(max_iters): d_next = d @ P # np.abs 取绝对值；.max() 取最大值。这是"无穷范数"——最大的那个差值。 if np.abs(d_next - d).max() < tol: return d_next, True d = d_next return d, False # ============================================================================= # 第 4 部分：从数据中学习转移矩阵 # ============================================================================= def estimate_transition_matrix(trajectories, n_states): """ 只给轨迹数据，通过"数数"估计转移矩阵（最大似然估计）。 【思路】 P_hat[s, s'] = （从 s 转移到 s' 的次数）/（从 s 出发的总次数） 这是讲义 7.7 节的内容，也是"基于模型的 RL"中"学习模型"这一步最简单的形式。 第 3 讲的"确定性等价"方法就是：先用这个方法估计 P（和 R），再把估计值当真，去做规划。 【参数】 trajectories : list[list[int]]，若干条轨迹，每条是一个状态列表。 例：[[3, 4, 5], [3, 2, 2, 1]] n_states : int，状态总数 N（不能从数据中推断，因为有的状态可能从没出现过）。 【返回】 (P_hat, counts) 二元组： - P_hat : 形状 (N, N) 的估计转移矩阵 - counts : 形状 (N, N) 的计数矩阵，counts[s, s'] = 观察到 s→s' 的次数 【边界】（重要！） - **从未出发过的状态**：分母为 0，无法估计。这里的处理是给它一个均匀分布 （"一无所知时，假设每种可能性一样"），保证 P_hat 仍是合法的转移矩阵。 但这只是一个"占位"，调用者应通过 counts.sum(axis=1) 检查每行的样本量， 样本量为 0 或很少的行都是不可信的。 这个问题在第 10 讲（离线 RL：数据没覆盖的状态怎么办）和 第 13 讲（探索：要不要主动去没去过的状态）中是核心问题。 - 轨迹中出现 >= n_states 或 < 0 的状态 → ValueError。 - 长度为 1 的轨迹没有任何转移，不贡献计数（合法，但没用）。 - 空列表 trajectories=[] → 所有行都是均匀分布。 """ counts = np.zeros((n_states, n_states)) for traj in trajectories: # 【Python 小知识】zip(traj[:-1], traj[1:]) 把"相邻两个元素"配成对。 # traj = [3, 4, 5, 4] # traj[:-1] = [3, 4, 5] （去掉最后一个） # traj[1:] = [4, 5, 4] （去掉第一个） # zip 后得到 (3,4), (4,5), (5,4) —— 正好是每一次转移 (s, s')。 for s, s_next in zip(traj[:-1], traj[1:]): if not (0 <= s < n_states and 0 <= s_next < n_states): raise ValueError(f"轨迹中出现了非法状态：({s}, {s_next})，状态必须在 0 到 {n_states - 1} 之间") counts[s, s_next] += 1 # keepdims=True 让结果保持形状 (N, 1) 而不是 (N,)， # 这样 counts / row_totals 时，每一行会自动除以该行自己的总数（这叫"广播 broadcasting"）。 row_totals = counts.sum(axis=1, keepdims=True) P_hat = np.empty((n_states, n_states)) # np.empty 创建矩阵但不初始化数值（马上会全部填上） visited = row_totals[:, 0] > 0 # 布尔向量：每个状态是否至少出发过一次 # 【Python 小知识】布尔索引：P_hat[visited] 只选出 visited 为 True 的那些行。 P_hat[visited] = counts[visited] / row_totals[visited] P_hat[~visited] = 1.0 / n_states # ~ 是"取反"：没访问过的行 → 均匀分布 return P_hat, counts # ============================================================================= # 第 5 部分：用数据验证马尔可夫性 # ============================================================================= def empirical_markov_check(trajectories, current_state, n_states): """ 用数据检验：给定"现在在 current_state"，"上一步在哪"会不会影响"下一步去哪"？ 【思路】 马尔可夫性说：P(s_{t+1} | s_t, s_{t-1}) = P(s_{t+1} | s_t)。 也就是说，只要知道现在在哪，过去在哪不重要。 我们对每一个可能的"上一步状态" prev，统计： 在 (s_{t-1}=prev, s_t=current_state) 的条件下，s_{t+1} 的经验分布。 如果这些分布都差不多（只有随机波动），说明马尔可夫性成立。 如果一个过程不是马尔可夫的，不同 prev 对应的分布会有**系统性**差异。 【参数】 trajectories : list[list[int]]，轨迹数据。 current_state : int，要检查的"现在"的状态。 n_states : int，状态总数。 【返回】 dict：{prev: (条件分布数组, 样本数)}，只包含在数据中出现过的 prev。 【边界】 - 某个 prev 样本很少时，它的经验分布波动会很大，不能据此下结论。 所以同时返回样本数，打印时一并展示。 - 轨迹长度 < 3 时没有"上一步-现在-下一步"三元组，不贡献数据。 """ counts = {} # prev -> 长度为 n_states 的计数数组 for traj in trajectories: # zip 三个错开的切片，得到连续的三元组 (s_{t-1}, s_t, s_{t+1}) for prev, cur, nxt in zip(traj[:-2], traj[1:-1], traj[2:]): if cur != current_state: continue # 【Python 小知识】continue 跳过本次循环剩下的代码，直接进入下一次 if prev not in counts: counts[prev] = np.zeros(n_states) counts[prev][nxt] += 1 result = {} # sorted(...) 把字典的键从小到大排序，让打印顺序固定、好读。 for prev in sorted(counts): total = counts[prev].sum() result[prev] = (counts[prev] / total, int(total)) return result # ============================================================================= # 第 6 部分：演示 # ============================================================================= def main(): rng = np.random.default_rng(seed=0) P = build_mars_rover_chain() # --------------------------------------------------------------------- print("=" * 78) print("演示 1：火星车马尔可夫链的转移矩阵（行 = 从哪来，列 = 到哪去）") print("=" * 78) print(" " + " ".join(f"s{j + 1}" for j in range(7))) for i in range(7): # 【Python 小知识】f"{x:.1f}" 保留 1 位小数；" ".join 把各列拼起来。 row = " ".join(f"{p:.1f}" for p in P[i]) print(f" s{i + 1}: {row} 行和 = {P[i].sum():.1f}") check_transition_matrix(P) print(" ✓ 通过合法性检查：方阵、非负、每行和为 1") # 故意构造一个错误的矩阵，看看检查函数如何报错 bad_P = P.copy() # .copy() 复制一份，避免修改原矩阵（直接 bad_P = P 只是起了个别名！） bad_P[2, 2] = 0.5 # 第 3 行和变成 0.4 + 0.5 + 0.4 = 1.3 try: check_transition_matrix(bad_P) except ValueError as e: print(f" ✗ 故意改错一个数后：ValueError: {e}") # --------------------------------------------------------------------- print("\n" + "=" * 78) print("演示 2：从 s4 出发采样 5 条轨迹（每条走 8 步）") print("=" * 78) for k in range(5): traj = sample_trajectory(P, start_state=3, n_steps=8, rng=rng) print(f" 轨迹 {k + 1}: {format_trajectory(traj)}") print(" → 同一个起点，每次走出的轨迹都不同：这就是\"随机过程\"。") # --------------------------------------------------------------------- print("\n" + "=" * 78) print("演示 3：状态分布的演化 d_t = d_0 P^t（从 s4 出发）") print("=" * 78) d0 = np.zeros(7) d0[3] = 1.0 # 确定在 s4 dists = propagate_distribution(P, d0, n_steps=50) print(" t " + "".join(f"{f's{j + 1}':>7s}" for j in range(7))) for t in [0, 1, 2, 5, 10, 50]: print(f" {t:2d} " + "".join(f"{x:7.3f}" for x in dists[t])) d_star, converged = stationary_distribution(P) print(f"\n 平稳分布 d* = {d_star} （收敛：{converged}）") print(f" 验证 d* P = d*：最大误差 = {np.abs(d_star @ P - d_star).max():.2e}") print(" → 无论从哪里出发，长期来看火星车在每个格子的概率都是 1/7 ≈ 0.143。") # 边界演示：周期链 P_periodic = np.array([[0.0, 1.0], [1.0, 0.0]]) d_osc = propagate_distribution(P_periodic, [1.0, 0.0], n_steps=4) print("\n 边界情况——周期链 P=[[0,1],[1,0]]，从 [1,0] 出发的分布序列：") for t, d in enumerate(d_osc): # 【Python 小知识】enumerate(列表) 同时给出下标和元素：(0, 第0个), (1, 第1个), ... print(f" t={t}: {d}") print(" → 永远在两个状态间来回振荡，分布不收敛。幂迭代并非对所有链都有效。") # --------------------------------------------------------------------- print("\n" + "=" * 78) print("演示 4：只看数据，把转移矩阵\"学\"出来（数数 = 最大似然估计）") print("=" * 78) for n_traj in [5, 50, 500, 5000]: data = [sample_trajectory(P, start_state=3, n_steps=20, rng=rng) for _ in range(n_traj)] P_hat, counts = estimate_transition_matrix(data, n_states=7) max_err = np.abs(P_hat - P).max() min_visits = int(counts.sum(axis=1).min()) print(f" {n_traj:5d} 条轨迹：最大估计误差 = {max_err:.3f}，访问最少的状态只出发过 {min_visits:6d} 次") print(" → 数据越多，估计越准。误差大约按 1/sqrt(样本数) 缩小：数据翻 100 倍，误差约缩小到 1/10。") print("\n 用 5000 条轨迹估计出的第 4 行（从 s4 出发），对比真实值：") print(f" 估计：{P_hat[3]}") print(f" 真实：{P[3]}") # 边界演示：有状态从未被访问 few = [[3, 4, 3]] # 只有一条很短的轨迹 P_few, counts_few = estimate_transition_matrix(few, n_states=7) print("\n 边界情况——只有一条轨迹 s4 → s5 → s4：") print(f" 各状态出发次数：{counts_few.sum(axis=1).astype(int)}") print(f" s1 那一行（从没去过）被填成均匀分布：{P_few[0]}") print(" → 这只是占位，完全不可信。\"没见过的状态怎么办\"是离线 RL 和探索的核心问题。") # --------------------------------------------------------------------- print("\n" + "=" * 78) print("演示 5：用数据验证马尔可夫性——上一步在哪，会影响下一步去哪吗？") print("=" * 78) data = [sample_trajectory(P, start_state=3, n_steps=50, rng=rng) for _ in range(2000)] result = empirical_markov_check(data, current_state=3, n_states=7) print(" 条件：现在在 s4。按\"上一步所在位置\"分组，统计下一步去向的经验分布：") print(" 去 s3 留 s4 去 s5 样本数") for prev, (dist, total) in result.items(): print(f" 上一步在 s{prev + 1}: {dist[2]:.3f} {dist[3]:.3f} {dist[4]:.3f} {total:6d}") print(f" 真实 P(·|s4): {P[3, 2]:.3f} {P[3, 3]:.3f} {P[3, 4]:.3f}") print(" → 三组分布几乎一样（差异只是随机波动）：\"上一步在哪\"不提供额外信息。") print(" 这正是马尔可夫性：给定现在，未来与过去无关。") if __name__ == "__main__": main()