""" ============================================================================= 第 1 讲 · 代码 3：马尔可夫奖励过程（MRP）——回报、折扣与价值的三种算法 ============================================================================= 对应讲义：lectures/lecture01_intro_to_rl.md 第 8 节 ----------------------------------------------------------------------------- 【这个文件要解决什么问题？】 ----------------------------------------------------------------------------- MRP = 马尔可夫链 + 奖励 + 折扣因子，即 (S, P, R, γ)。 本讲最重要的问题是：**每个状态的价值 V(s) 是多少？** V(s) = E[ r_t + γ r_{t+1} + γ² r_{t+2} + ... | s_t = s ] 我们用三种完全不同的思路来算它，并验证三者得到相同的答案： 方法一 蒙特卡洛 ：模拟很多条轨迹，把回报取平均。 （只需要能采样） 方法二 解析解 ：解线性方程组 (I - γP) V = R。 （需要知道 P 和 R） 方法三 迭代法 ：反复套用 V ← R + γ P V 直到不变。 （需要知道 P 和 R） ----------------------------------------------------------------------------- 【文件结构】 ----------------------------------------------------------------------------- 第 1 部分 构建火星车 MRP 第 2 部分 回报的计算（正着算、倒着算） 第 3 部分 方法一：蒙特卡洛 第 4 部分 方法二：解析解 第 5 部分 方法三：迭代法 第 6 部分 演示：复现讲义中的所有数字 ----------------------------------------------------------------------------- 【运行方式】 ----------------------------------------------------------------------------- python 03_markov_reward_process.py 【下标约定】代码中的状态 i <==> 讲义中的 s_{i+1} 【奖励约定】处于状态 s_t 时获得奖励 r_t = R(s_t)，然后转移到 s_{t+1}（CS234 约定） """ import numpy as np np.set_printoptions(precision=2, suppress=True) # ============================================================================= # 第 1 部分：构建火星车 MRP # ============================================================================= def build_mars_rover_mrp(): """ 构建讲义 8.1 节的火星车 MRP（不含 γ，γ 在调用各算法时单独传入）。 【返回】 (P, R) 二元组： - P : 形状 (7, 7) 的转移矩阵（与 02_markov_chain.py 相同） - R : 长度为 7 的奖励向量，R[s] = 在状态 s 获得的期望即时奖励 【为什么 γ 不放在这里？】 讲义 8.4 节说过：γ 是"问题定义"的一部分，不同的 γ 代表不同的"耐心程度"。 把它作为参数单独传入，方便我们用同一个 (P, R) 比较不同 γ 的效果。 """ n = 7 P = np.zeros((n, n)) P[0, 0], P[0, 1] = 0.6, 0.4 # 【Python 小知识】a, b = 1, 2 可以同时给两个变量赋值。 for i in range(1, n - 1): P[i, i - 1], P[i, i], P[i, i + 1] = 0.4, 0.2, 0.4 P[n - 1, n - 2], P[n - 1, n - 1] = 0.4, 0.6 R = np.array([1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 10.0]) return P, R def _check_gamma(gamma, allow_one=True): """ 检查折扣因子是否合法。（函数名以下划线开头，是 Python 的惯例，表示"内部辅助函数"。） 【参数】 gamma : float，折扣因子。 allow_one : bool，是否允许 γ = 1。 - 计算有限长度序列的回报时，γ = 1 完全没问题（就是不打折的求和）； - 解析解要求 γ < 1，否则矩阵 I - P 不可逆。 【边界】 γ < 0 或 γ > 1 都没有意义： - γ < 0 会让奖励的符号正负交替，毫无解释； - γ > 1 表示"越远的未来越重要"，无限视界时回报会爆炸。 """ upper_ok = gamma <= 1.0 if allow_one else gamma < 1.0 if not (0.0 <= gamma and upper_ok): valid = "[0, 1]" if allow_one else "[0, 1)" raise ValueError(f"gamma 必须在 {valid} 内，但收到了 {gamma}") # ============================================================================= # 第 2 部分：回报（Return）的计算 # ============================================================================= def discounted_return(rewards, gamma): """ 计算一条奖励序列的折扣回报（"正着算"）： G_0 = r_0 + γ r_1 + γ² r_2 + ... + γ^{T-1} r_{T-1} 【思路】 最直接的写法：第 k 个奖励乘以 γ^k，然后全部加起来。 【参数】 rewards : 一串数字（list 或 numpy 数组），[r_0, r_1, ..., r_{T-1}]。 gamma : float，折扣因子，范围 [0, 1]。 【返回】 float，回报 G_0。 【边界】 - 空序列 [] → 返回 0.0（没有奖励，回报自然是 0）。 - γ = 0 → 只剩 r_0（注意 Python 中 0 ** 0 = 1，所以第 0 项的权重仍是 1，正确）。 - γ = 1 → 普通求和。 【例子】 discounted_return([0, 0, 0, 10], 0.5) → 0 + 0 + 0 + 0.125 × 10 = 1.25 """ _check_gamma(gamma, allow_one=True) rewards = np.asarray(rewards, dtype=float) # 【Python 小知识】 # np.arange(T) 生成 [0, 1, ..., T-1]； # gamma ** np.arange(T) 对每个元素求幂，得到 [1, γ, γ², ..., γ^{T-1}]（这叫"向量化"运算， # 不用写 for 循环，numpy 会一次性对整个数组计算，又快又简洁）； # np.sum(a * b) 是两个向量逐元素相乘再求和（即点积）。 discounts = gamma ** np.arange(len(rewards)) return float(np.sum(discounts * rewards)) def returns_for_every_timestep(rewards, gamma): """ 一次性算出序列中**每一个时刻**的回报 G_0, G_1, ..., G_{T-1}（"倒着算"）。 【思路】 如果对每个 t 都用 discounted_return 重新求和，总共要 O(T²) 次运算。 但讲义 8.3 节的递归关系告诉我们： G_t = r_t + γ · G_{t+1}， 且 G_T = 0（序列结束后没有奖励） 所以可以**从最后一个时刻往前**扫一遍，每步只做一次乘法和一次加法，总共 O(T)。 以 rewards = [0, 0, 10, 10, 0], γ = 0.9 为例（讲义练习 3）： G_4 = 0 + 0.9 × 0 = 0 G_3 = 10 + 0.9 × 0 = 10 G_2 = 10 + 0.9 × 10 = 19 G_1 = 0 + 0.9 × 19 = 17.1 G_0 = 0 + 0.9 × 17.1 = 15.39 【RL 联系】 这个"倒着算回报"的技巧在第 3 讲蒙特卡洛方法、第 6 讲 REINFORCE 中都会直接用到。 而 "G_t = r_t + γ G_{t+1}" 这个关系本身，就是贝尔曼方程和 TD 学习的源头。 【参数】 rewards : 一串数字 [r_0, ..., r_{T-1}]。 gamma : float，范围 [0, 1]。 【返回】 numpy 数组，长度 T，第 t 个元素是 G_t。 【边界】 - 空序列 → 返回空数组。 """ _check_gamma(gamma, allow_one=True) rewards = np.asarray(rewards, dtype=float) T = len(rewards) # np.zeros(T) 预先分配好结果数组 returns = np.zeros(T) G = 0.0 # 对应 G_T = 0：最后一步之后没有奖励 # 【Python 小知识】range(T - 1, -1, -1) 生成 T-1, T-2, ..., 0（倒序）。 # range(起点, 终点, 步长)，终点不包含，所以终点写 -1 才能包含 0。 for t in range(T - 1, -1, -1): G = rewards[t] + gamma * G returns[t] = G return returns # ============================================================================= # 第 3 部分：方法一 —— 蒙特卡洛 # ============================================================================= def recommended_horizon(gamma, r_max, tolerance): """ 蒙特卡洛需要把无限长的轨迹截断到有限步 H。这个函数算出 H 至少取多大， 才能保证截断带来的误差不超过 tolerance。 【思路】 截断后丢掉的是第 H 步及以后的所有奖励，它们的折扣和最多为： γ^H · r_max + γ^{H+1} · r_max + ... = γ^H · r_max / (1 - γ) 要让它 <= tolerance，解出： γ^H <= tolerance · (1 - γ) / r_max H >= log(tolerance · (1 - γ) / r_max) / log(γ) （注意 log(γ) < 0，两边除以负数时不等号方向要反过来，所以是 >=。） 【参数】 gamma : float，范围 [0, 1)。 r_max : float，奖励绝对值的上界，必须 > 0。 tolerance : float，允许的截断误差，必须 > 0。 【返回】 int，推荐的视界 H（至少为 1）。 【边界】 - γ = 0：只有第一步奖励有用，H = 1 就够。 - γ = 1：丢掉的尾巴没有折扣，无论 H 多大截断误差都可能不小 → ValueError。 - 如果 tolerance 很大，算出来的 H 可能 <= 0，我们至少返回 1。 """ _check_gamma(gamma, allow_one=False) if r_max <= 0 or tolerance <= 0: raise ValueError("r_max 和 tolerance 都必须 > 0") if gamma == 0.0: return 1 # np.ceil 向上取整：比如算出 H >= 12.3，就取 13 H = np.log(tolerance * (1 - gamma) / r_max) / np.log(gamma) return max(1, int(np.ceil(H))) def sample_episode_rewards(P, R, start_state, horizon, rng): """ 从 start_state 出发，采样一条长度为 horizon 的轨迹，返回途经的状态和获得的奖励。 【参数】 P, R : MRP 的转移矩阵和奖励向量。 start_state : int，起点。 horizon : int，步数 H，必须 >= 1。 rng : numpy 随机数生成器。 【返回】 (states, rewards)： - states : list[int]，长度 H，[s_0, ..., s_{H-1}] - rewards : list[float]，长度 H，[r_0, ..., r_{H-1}]，其中 r_t = R[s_t] 【注意】 这里我们只记录 H 个状态（在 s_{H-1} 拿完奖励就停）， 所以 H 步只需要 H-1 次转移。和 02 文件中 sample_trajectory 的"走 n 步经过 n+1 个状态" 是同一件事的两种记法，写代码时要时刻清楚自己用的是哪种。 """ if horizon < 1: raise ValueError(f"horizon 必须 >= 1，但收到了 {horizon}") n = P.shape[0] states, rewards = [], [] s = start_state for t in range(horizon): states.append(s) rewards.append(float(R[s])) if t < horizon - 1: # 最后一步拿完奖励就结束，不需要再转移 s = int(rng.choice(n, p=P[s])) return states, rewards def monte_carlo_value(P, R, gamma, n_episodes, horizon, rng): """ 【方法一】用蒙特卡洛模拟估计每个状态的价值。 【思路】 V(s) 的定义就是"从 s 出发的回报的期望"。期望算不出来？那就采样！ 对每个状态 s： 1. 从 s 出发模拟 n_episodes 条轨迹，每条长 horizon 步； 2. 算出每条轨迹的折扣回报 G； 3. 把这些 G 取平均，作为 V(s) 的估计。 依据是**大数定律**：样本平均值会收敛到期望值。 【RL 联系】 注意这个函数里 P 只被用来"采样下一个状态"（通过 sample_episode_rewards）。 如果把采样换成"与真实环境交互"，我们就完全不需要知道 P 的数值！ 这正是第 3 讲"无模型策略评估"中蒙特卡洛方法的核心思想。 【参数】 P, R : MRP。 gamma : float，范围 [0, 1]。 n_episodes : int，每个状态模拟的轨迹数，必须 >= 1。越大越准，误差 ∝ 1/sqrt(n_episodes)。 horizon : int，每条轨迹的长度 H。应当用 recommended_horizon 选取， 否则截断误差会让估计系统性地偏小（奖励为正时）。 rng : numpy 随机数生成器。 【返回】 (V_hat, std_errors)： - V_hat : 长度 N 的数组，价值估计 - std_errors : 长度 N 的数组，每个估计的标准误差（衡量随机误差的大小） 【边界】 - γ = 1 时截断误差不会随 H 衰减，估计的其实是"H 步有限视界"的价值，而不是无限视界价值。 - n_episodes = 1 时无法估计标准误差，返回 0。 - 计算量 = N × n_episodes × horizon，状态多、γ 接近 1 时会很慢。 """ _check_gamma(gamma, allow_one=True) if n_episodes < 1: raise ValueError(f"n_episodes 必须 >= 1，但收到了 {n_episodes}") n = P.shape[0] V_hat = np.zeros(n) std_errors = np.zeros(n) for s in range(n): returns = [] for _ in range(n_episodes): _, rewards = sample_episode_rewards(P, R, s, horizon, rng) returns.append(discounted_return(rewards, gamma)) returns = np.array(returns) V_hat[s] = returns.mean() if n_episodes > 1: std_errors[s] = returns.std(ddof=1) / np.sqrt(n_episodes) return V_hat, std_errors # ============================================================================= # 第 4 部分：方法二 —— 解析解 # ============================================================================= def analytic_value(P, R, gamma): """ 【方法二】直接解贝尔曼方程的矩阵形式，得到精确的价值。 【思路】 MRP 贝尔曼方程（讲义 8.7 节推导）： V(s) = R(s) + γ Σ_{s'} P(s'|s) V(s') 对所有状态写在一起，就是矩阵形式： V = R + γ P V 移项： (I - γP) V = R 这是一个标准的线性方程组 "A x = b"，其中 A = I - γP，x = V，b = R。 【为什么用 np.linalg.solve 而不是 np.linalg.inv？】 数学上 V = (I - γP)^{-1} R，但在数值计算中， "先求逆矩阵再相乘"比"直接解方程组"更慢、误差更大。 经验法则：**只要你的目的是解方程，就用 solve，永远不要显式求逆。** 【参数】 P : 转移矩阵，形状 (N, N)。 R : 奖励向量，长度 N。 gamma : float，范围 [0, 1)。**注意不能等于 1。** 【返回】 长度 N 的数组，精确的价值 V。 【边界】 - γ = 1 → ValueError。原因：P 的每行和为 1，所以 (I - P) 乘以全 1 向量等于 0 向量， 矩阵 I - P 奇异（不可逆），方程要么无解要么有无穷多解。 直觉上：不打折又永不终止，价值可能是无穷大。 （如果 MRP 有"终止状态"——到达后不再有奖励且不再转移——γ = 1 也可以求解， 但需要把终止状态单独处理，本讲不涉及。） - 计算量约 O(N³)：N = 1000 时很快，N = 100 万时就不可行了。 - P 和 R 的尺寸不匹配 → ValueError。 """ _check_gamma(gamma, allow_one=False) n = P.shape[0] R = np.asarray(R, dtype=float) if P.shape != (n, n) or R.shape != (n,): raise ValueError(f"尺寸不匹配：P 的形状 {P.shape}，R 的形状 {R.shape}") # np.eye(n) 是 n×n 单位矩阵（对角线为 1，其余为 0）。 A = np.eye(n) - gamma * P return np.linalg.solve(A, R) # ============================================================================= # 第 5 部分：方法三 —— 迭代法（动态规划） # ============================================================================= def iterative_value(P, R, gamma, tol=1e-10, max_iters=100_000): """ 【方法三】反复套用贝尔曼方程，让价值估计逐步逼近真实值。 【思路】 从 V_0 = 0 开始，不断执行： V_k = R + γ P V_{k-1} 直到相邻两次的变化足够小。 V_k 的直觉含义："如果只看未来 k 步，能拿多少奖励"： V_1 = R （只看 1 步） V_2 = R + γ P R （看 2 步） ... k 越大看得越远，远处的奖励被 γ^k 压得越来越小，所以 V_k 收敛。 【为什么一定收敛？（讲义 8.7 节的压缩映射）】 设真实价值为 V，则 V_k - V = γ P (V_{k-1} - V)。 P 每一行都是加权平均，不会放大误差的最大值；再乘 γ < 1，误差至少缩小到 γ 倍。 本函数返回每次迭代的变化量 deltas，你可以亲眼看到 deltas[k+1]/deltas[k] ≈ γ。 【参数】 P, R : MRP。 gamma : float，范围 [0, 1]。允许传入 1，但通常不会收敛（见边界）。 tol : float，停止阈值。当 max_s |V_k(s) - V_{k-1}(s)| < tol 时停止。 max_iters : int，最多迭代次数，防止不收敛时死循环。 【返回】 (V, n_iters, deltas, converged)： - V : 长度 N 的数组，最后一次迭代的价值 - n_iters : int，实际迭代次数 - deltas : list[float]，每次迭代的最大变化量（用来观察收敛速度） - converged : bool，是否在 max_iters 次内收敛 【边界】 - γ = 1：对火星车这种永不终止、奖励为正的链，V_k 会无限增长，永远不收敛， 函数会在 max_iters 次后返回 converged=False。调用者**必须检查 converged**。 - γ 越接近 1 收敛越慢：迭代次数大约 ∝ 1/(1-γ)。 - 小心"收敛"的含义：变化量 < tol 不代表与真实值的误差 < tol。 可以证明真实误差 <= γ/(1-γ) × 变化量，γ 接近 1 时这个倍数很大。 """ _check_gamma(gamma, allow_one=True) n = P.shape[0] V = np.zeros(n) deltas = [] for k in range(1, max_iters + 1): # 这一行就是整个算法！R + γ P V，其中 P @ V 是矩阵乘向量： # (P @ V)[s] = Σ_{s'} P[s, s'] V[s'] = "从 s 出发，下一步所在状态的平均价值" V_new = R + gamma * (P @ V) delta = float(np.abs(V_new - V).max()) deltas.append(delta) V = V_new if delta < tol: return V, k, deltas, True return V, max_iters, deltas, False # ============================================================================= # 第 6 部分：演示 # ============================================================================= def pad_to_width(text, width): """ 把字符串在右侧补空格，使它在终端里显示的宽度恰好为 width。 【为什么需要它？】 终端里一个中文字符占 2 个英文字符的宽度，但 Python 的 len("蒙特卡洛") 是 4， 所以 f"{text:<14s}" 这种按字符数对齐的写法，遇到中文就会歪。 这里用 unicodedata.east_asian_width 判断每个字符是"宽字符"（W/F，如中文）还是普通字符。 【边界】 如果 text 本身已经比 width 宽，不截断，原样返回。 """ import unicodedata # 标准库，无需安装；放在函数内导入，表示只有这里用到它 display = sum(2 if unicodedata.east_asian_width(ch) in ("W", "F") else 1 for ch in text) return text + " " * max(0, width - display) def print_value_row(label, V): """把一个价值向量按讲义表格的格式打印成一行（label 占 16 列宽，每个数占 9 列）。""" print(" " + pad_to_width(label, 16) + "".join(f"{v:9.2f}" for v in V)) def main(): P, R = build_mars_rover_mrp() rng = np.random.default_rng(seed=0) header = " " + " " * 16 + "".join(f"{f's{j + 1}':>9s}" for j in range(7)) # --------------------------------------------------------------------- print("=" * 78) print("演示 1：手算回报（复现讲义 8.6 节，γ = 0.5，从 s4 出发的 3 条轨迹）") print("=" * 78) examples = [ ("s4 → s5 → s6 → s7", [3, 4, 5, 6]), ("s4 → s4 → s5 → s4", [3, 3, 4, 3]), ("s4 → s3 → s2 → s1", [3, 2, 1, 0]), ] for desc, states in examples: # 【Python 小知识】[float(R[s]) for s in states] 取出每个状态的奖励，组成奖励序列。 # 套一层 float(...) 是为了打印好看：numpy 2.x 直接打印 R[s] 会显示成 np.float64(0.0)。 rewards = [float(R[s]) for s in states] G = discounted_return(rewards, gamma=0.5) print(f" {desc} 奖励 {rewards} G_0 = {G}") # --------------------------------------------------------------------- print("\n" + "=" * 78) print("演示 2：倒着算每个时刻的回报（讲义练习 3：奖励 [0,0,10,10,0]，γ = 0.9）") print("=" * 78) rewards = [0, 0, 10, 10, 0] Gs = returns_for_every_timestep(rewards, gamma=0.9) for t, G in enumerate(Gs): print(f" G_{t} = {G:.2f}") print(f" 检查：正着算的 G_0 = {discounted_return(rewards, 0.9):.2f}，与倒着算一致 ✓") # --------------------------------------------------------------------- print("\n" + "=" * 78) print("演示 3：三种方法计算 V（γ = 0.5）") print("=" * 78) gamma = 0.5 H = recommended_horizon(gamma, r_max=10.0, tolerance=1e-3) print(f" 蒙特卡洛的截断视界：为使截断误差 <= 0.001，取 H = {H} 步\n") V_mc, se = monte_carlo_value(P, R, gamma, n_episodes=5000, horizon=H, rng=rng) V_exact = analytic_value(P, R, gamma) V_iter, n_iters, deltas, converged = iterative_value(P, R, gamma) print(header) print_value_row("蒙特卡洛", V_mc) print_value_row(" ±2×标准误差", 2 * se) print_value_row("解析解", V_exact) print_value_row("迭代法", V_iter) print(f"\n 迭代法用了 {n_iters} 次迭代收敛（converged = {converged}）") print(f" 解析解 vs 迭代法 最大差异：{np.abs(V_exact - V_iter).max():.2e}") print(f" 解析解 vs 蒙特卡洛 最大差异：{np.abs(V_exact - V_mc).max():.3f}（随机误差，样本越多越小）") print("\n 观察：V(s1)=1.53 < V(s6)=3.59 —— s6 自己没有奖励，但它紧挨着 +10，\"前途\"更好。") print(" 这就是讲义里说的\"价值 ≠ 奖励\"。") # --------------------------------------------------------------------- print("\n" + "=" * 78) print("演示 4：迭代法的收敛过程——每次迭代，变化量缩小到约 γ 倍") print("=" * 78) for g in [0.5, 0.9, 0.99]: _, n_it, ds, ok = iterative_value(P, R, g) # 取最后两次迭代的变化量之比（前几次迭代还没进入"稳定阶段"，比值会偏离 γ） # 【Python 小知识】ds[-1] 是列表最后一个元素，ds[-2] 是倒数第二个。 ratio = ds[-1] / ds[-2] print(f" γ = {g:<5} 迭代 {n_it:5d} 次收敛；最后两次变化量之比 = {ratio:.4f}（≈ γ）") print(" → γ 越接近 1，每次误差缩小得越少，需要的迭代次数越多（约 ∝ 1/(1-γ)）。") print("\n γ = 0.5 时前 8 次迭代的变化量：") for k in range(8): print(f" 第 {k + 1} 次：{deltas[k]:.6f}") # --------------------------------------------------------------------- print("\n" + "=" * 78) print("演示 5：γ 如何改变价值（复现讲义 8.8 节表格）") print("=" * 78) print(header + " 有效视界≈1/(1-γ)") for g in [0.0, 0.5, 0.9, 0.99]: V = analytic_value(P, R, g) print(" " + pad_to_width(f"γ = {g}", 16) + "".join(f"{v:9.2f}" for v in V) + f" {1 / (1 - g):8.0f} 步") print("\n 观察：") print(" - γ = 0：价值 = 即时奖励（完全短视）。") print(" - γ 变大：离 s7 近的状态价值上升更快，看得越远，+10 的\"引力\"越大。") print(" - γ = 0.99：价值很大，但各状态间的差距相对变小（长期看从哪出发区别不大）。") # --------------------------------------------------------------------- print("\n" + "=" * 78) print("演示 6：边界情况") print("=" * 78) try: analytic_value(P, R, gamma=1.0) except ValueError as e: print(f" [解析解 γ=1] ValueError: {e}") print(" （I - P 不可逆：不打折又永不终止，价值是无穷大）") V1, n1, ds1, ok1 = iterative_value(P, R, gamma=1.0, max_iters=200) print(f" [迭代法 γ=1] 迭代 {n1} 次后 converged = {ok1}，V(s4) 已经涨到 {V1[3]:.1f}，还在继续涨") try: discounted_return([1, 2, 3], gamma=1.5) except ValueError as e: print(f" [γ 越界] ValueError: {e}") print(f" [空序列] discounted_return([], 0.9) = {discounted_return([], 0.9)}") print(f" [γ=0] discounted_return([5, 100, 100], 0.0) = {discounted_return([5, 100, 100], 0.0)}（只看第一步）") # 截断视界太短会怎样？ V_short, _ = monte_carlo_value(P, R, 0.9, n_episodes=2000, horizon=5, rng=rng) V_true = analytic_value(P, R, 0.9) print(f" [视界太短] γ=0.9、H=5 时蒙特卡洛估计 V(s7) = {V_short[6]:.2f}，真实值 = {V_true[6]:.2f}") print(f" （截断丢掉了大量未来奖励，系统性偏小。推荐 H = " f"{recommended_horizon(0.9, 10.0, 1e-3)}）") if __name__ == "__main__": main()