{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 从马尔可夫过程（MP）到马尔可夫奖励过程（MRP），再到 MDP 的占用度量（occupancy measure）。\n",
    "\n",
    "1) MP：三状态天气链，求平稳分布；2) MRP：加上奖励，用 v = (I - γP)^{-1} r 求价值并用蒙特卡洛验证；\n",
    "3) 占用度量：在本章的两状态 MDP（A/B）上计算 ρ^π(s, a)，并验证 V^π(s0) = <ρ^π, r> / (1 - γ)。\n",
    "\n",
    "配套教程：[马尔可夫决策过程（MDP）](https://dqn.cn/basics/mdp)（DQN.CN 强化学习中文教程）\n",
    "\n",
    "源码：`examples/basics/mdp_occupancy.py`，本 notebook 由 `scripts/build-notebooks.mjs` 自动生成，代码与 .py 文件一致。代码默认在 CPU 上运行，固定了随机种子；部分脚本支持设置环境变量 `DQN_QUICK=1` 只跑很短的版本（见下一个单元）。"
   ],
   "id": "basics-mdp_occupancy-py-0"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 环境准备：在 Colab 等缺少依赖的环境中自动安装（本地已装好依赖时什么也不做）\n",
    "import importlib.util, os, subprocess, sys\n",
    "sys.argv = sys.argv[:1]   # notebook 内核会带上 -f 等参数，清掉以免脚本把它们当成命令行参数\n",
    "# os.environ[\"DQN_QUICK\"] = \"1\"   # 取消注释：支持该开关的脚本只跑很短的版本，用于快速检查"
   ],
   "id": "basics-mdp_occupancy-py-1"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "rng = np.random.default_rng(0)\n",
    "gamma = 0.9\n",
    "\n",
    "# ---------- 1) 马尔可夫过程：状态 0=晴 1=阴 2=雨 ----------\n",
    "names = [\"晴\", \"阴\", \"雨\"]\n",
    "P = np.array([[0.7, 0.2, 0.1],\n",
    "              [0.3, 0.4, 0.3],\n",
    "              [0.2, 0.4, 0.4]])\n",
    "w, V = np.linalg.eig(P.T)                                    # 平稳分布：P^T μ = μ 的特征向量\n",
    "mu = np.real(V[:, np.argmin(np.abs(w - 1))]); mu /= mu.sum()\n",
    "print(\"MP 平稳分布：\", {n: round(float(m), 4) for n, m in zip(names, mu)})\n",
    "s, counts = 0, np.zeros(3)\n",
    "for _ in range(200_000):\n",
    "    s = rng.choice(3, p=P[s]); counts[s] += 1\n",
    "print(\"MP 模拟 20 万步的状态频率：\", {n: round(float(c), 4) for n, c in zip(names, counts / counts.sum())})\n",
    "\n",
    "# ---------- 2) 马尔可夫奖励过程：每天的“出行满意度”奖励 ----------\n",
    "r = np.array([2.0, 0.5, -1.0])\n",
    "v = np.linalg.solve(np.eye(3) - gamma * P, r)\n",
    "print(\"MRP 解析价值 v = (I-γP)^{-1} r：\", {n: round(float(x), 3) for n, x in zip(names, v)})\n",
    "mc = []\n",
    "for _ in range(5000):                                        # 从“雨”出发的蒙特卡洛估计（截断 150 步）\n",
    "    s, G = 2, 0.0\n",
    "    for t in range(150):\n",
    "        G += gamma ** t * r[s]; s = rng.choice(3, p=P[s])\n",
    "    mc.append(G)\n",
    "print(f\"MRP 蒙特卡洛估计 v(雨)：{np.mean(mc):.3f} ± {np.std(mc) / np.sqrt(len(mc)):.3f}\")\n",
    "\n",
    "# ---------- 3) 占用度量：两状态 MDP，状态 A=0, B=1；A 的动作 0=停留 1=前往，B 只有一个动作 ----------\n",
    "# 为了写成统一的矩阵形式，B 的两个动作效果相同（奖励 2，留在 B）\n",
    "Pmdp = np.zeros((2, 2, 2)); R = np.zeros((2, 2))\n",
    "Pmdp[0, 0, 0], R[0, 0] = 1.0, 1.0                            # A 停留：奖励 1，留在 A\n",
    "Pmdp[0, 1, 1], R[0, 1] = 1.0, 0.0                            # A 前往：奖励 0，到 B\n",
    "Pmdp[1, :, 1], R[1, :] = 1.0, 2.0                            # B：奖励 2，留在 B\n",
    "s0 = np.array([1.0, 0.0])                                    # 初始总在 A"
   ],
   "id": "basics-mdp_occupancy-py-2"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def occupancy(pi):\n",
    "    \"\"\"ρ(s,a) = (1-γ) Σ_t γ^t Pr(s_t=s, a_t=a)。先解 d = (1-γ) s0 + γ P_π^T d，再乘 π。\"\"\"\n",
    "    P_pi = np.einsum(\"sa,sat->st\", pi, Pmdp)\n",
    "    d = np.linalg.solve(np.eye(2) - gamma * P_pi.T, (1 - gamma) * s0)\n",
    "    return d[:, None] * pi"
   ],
   "id": "basics-mdp_occupancy-py-3"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def value(pi):\n",
    "    P_pi = np.einsum(\"sa,sat->st\", pi, Pmdp); r_pi = (pi * R).sum(1)\n",
    "    return np.linalg.solve(np.eye(2) - gamma * P_pi, r_pi)"
   ],
   "id": "basics-mdp_occupancy-py-4"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for name, pA in [(\"总是停留\", 0.0), (\"总是前往\", 1.0), (\"随机 50/50\", 0.5)]:\n",
    "    pi = np.array([[1 - pA, pA], [0.5, 0.5]])\n",
    "    rho = occupancy(pi)\n",
    "    print(f\"策略「{name}」：ρ(A,停留)={rho[0, 0]:.3f} ρ(A,前往)={rho[0, 1]:.3f} ρ(B,·)={rho[1].sum():.3f} \"\n",
    "          f\"| Σρ={rho.sum():.3f} | <ρ,r>/(1-γ)={(rho * R).sum() / (1 - gamma):.3f} | V(A)={value(pi)[0]:.3f}\")"
   ],
   "id": "basics-mdp_occupancy-py-5"
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python"
  },
  "colab": {
   "provenance": []
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}