# 教材位置：第10章，第10.5节；图 10-3。
# 单位、变量与形状：P 为 (3,3)，x0 为 (3,)，pred 为 (2000,3)；面积占比无量纲，时间步按年
# 计。
# 案例 35：土地类型马尔可夫转移、期望占比与转移矩阵的情景不确定性。
# 学习目标：区分等面积单元的多项抽样、固定矩阵的十步期望，以及扰动转移矩阵后的情景分布。
# 核心关系：行向量 xₜ₊₁=xₜP，固定矩阵十步期望为 x₀P¹⁰；稳态满足 πP=π。
# 数据与形状：三类土地，P 为 3×3，x0 为长度 3 的初始面积比例；
# 各行分别抽 狄利克雷(300*P[i]) 得到情景矩阵 Q，每个情景内固定 Q 连续推演 10 步。
# 2000 份情景结果 pred 为 (2000,3)，报告各分量的 5%、50%、95% 边际分位数。
# 这些是给定扰动设定的情景分位区间，不是用观测样本推导的频率学置信区间；
# 也未指定后验数据模型，因此不能自动称为后验可信区间或真实土地变化预测区间。
# 输出：抽样计数、稳态、十步点值和情景分位数；生成 application_land.png。
# 运行方式：安装依赖后在本目录执行 python e35_land.py，图片保存到当前工作目录。

# NumPy 用于多项/狄利克雷 抽样、矩阵幂、特征分解以及情景分位数。
import numpy as np
# 左图比较点预测和情景区间，右图展示每类土地的边际情景直方图。
import matplotlib.pyplot as plt
# 在创建第一幅图之前配置可跨系统选择的中文字体。
from plot_style import configure_chinese_plot

configure_chinese_plot()
# 固定种子 2026，先抽一次多项计数，再依次抽取 2000 个转移矩阵情景。
rng = np.random.default_rng(2026)
# 行与列的类型顺序均为农地、建设用地、其他用地；P[i,j] 是从类型 i 到 j 的一步概率。
# 所有元素都为正，使本例固定矩阵具有唯一稳态且长期期望分布收敛。
P = np.array([[0.88, 0.09, 0.03], [0.04, 0.93, 0.03], [0.02, 0.05, 0.93]])
# 初始三类土地面积占比分别为 42%、33%、25%；作为行概率向量右乘转移矩阵。
x0 = np.array([0.42, 0.33, 0.25])
# 验证每行是有效概率向量，保证乘法更新保持总面积比例为 1。
assert np.all(P >= 0) and np.allclose(P.sum(axis=1), 1)
# 初始占比也必须非负且和为 1，否则不能解释为三类土地的完整组成。
assert np.all(x0 >= 0) and np.isclose(x0.sum(), 1)
# 单独示范：20000 个当前都属于第 0 类的等面积单元，下一步去向按 P[0] 多项抽样。
# 它不是按 x0 配置的全区域单元样本，也没有被用来拟合后面的 狄利克雷 参数。
cell_counts = rng.multinomial(20000, P[0])
# 打印三种去向的实现计数；其期望为 20000*P[0]，本次抽样一般有随机波动。
print("等面积单元下一年的去向计数：", cell_counts)
# 多项抽样把所有单元分到三个互斥且完备的去向，因此计数总和必须恰为 20000。
assert cell_counts.sum() == 20000
# πP=π 转置成 Pᵀπᵀ=πᵀ，所以用 P.T 的右特征向量计算稳态。
vals, vecs = np.linalg.eig(P.T)
# 选择最接近特征值 1 的那一列；本例理论稳态为实数，取实部去除数值虚部。
# 特征向量尺度任意，除以分量和后才能解释为概率占比；steady 只对应固定的 P。
steady = np.real(vecs[:, np.argmin(np.abs(vals - 1))])
steady /= steady.sum()
# 检查归一化结果非负且满足固定矩阵的稳态方程，防止特征向量选错或方向用错。
assert np.all(steady >= 0) and np.allclose(steady @ P, steady)
# 计算固定 P 下第 10 步的期望面积占比，不是抽样得到的实际面积计数。
# 有限步点值与长期稳态概念不同，十步时不必已经达到 steady。
deterministic = x0 @ np.linalg.matrix_power(P, 10)
# 先用列表逐份收集情景结果，循环结束再转换为便于按列统计的二维数组。
pred = []
# 构造 2000 个独立的转移矩阵情景；下划线说明无需使用循环编号。
for _ in range(2000):
    # 狄利克雷 每次产生一行和为 1 的非负向量；vstack 把三行拼成 3×3 行随机矩阵。
    # 各行参数总和为 300，均值正好等于相应 P 行，单个分量方差为 pⱼ(1-pⱼ)/301。
    # 300 是人为设定的集中程度，越大扰动通常越小；没有在本例中从真实转移计数估计它。
    # 同一行的分量因总和固定而相关，不能视为三个独立概率。
    Q = np.vstack([rng.dirichlet(300 * row) for row in P])
    # 每份情景先抽一个 Q，再把这个 Q 固定使用 10 步，表达时间齐次转移参数的情景变化。
    # 此处没有逐年重新抽 Q，也没有模拟有限单元数的逐步状态抽样，结果仍是给定 Q 的期望占
    # 比。
    pred.append(x0 @ np.linalg.matrix_power(Q, 10))
# 转换为 (2000,3) 数组：行是情景，列是第 10 步的三种土地类型占比。
pred = np.asarray(pred)
# 验证每一份情景结果都在概率单纯形上，即非负且三类比例合计为 1。
assert np.all(pred >= 0) and np.allclose(pred.sum(axis=1), 1)
# 沿情景行求每一类的 5%、中位数、95% 分位数，三个输出各为长度 3 的数组。
# 每一类的 lo～hi 是中央约 90% 的边际情景区间，不是均值置信区间。
# 三类边际中位数一般不保证相加为 1，各列分位端点也不代表同一份可共同实现的情景，
# 因此不能把三条边际区间解释成联合覆盖率 90% 的整体区域。
lo, mid, hi = np.quantile(pred, [0.05, 0.5, 0.95], axis=0)
# 报告固定 P 的长期稳态比例，显示到 4 位小数。
print("固定转移矩阵的稳态占比：", np.round(steady, 4))
# 报告固定 P 的有限十步点值，便于与稳态和情景中位数分别比较。
print("固定转移矩阵的十步期望占比：", np.round(deterministic, 4))
# 按三行输出 5%、50%、95% 分位数；每列对应一种土地类型。
print("边际情景分位数（5%、50%、95%）：\n",
      np.round([lo, mid, hi], 4))
# 创建情景分位区间图与情景分布图，两图使用同一份 pred 结果。
fig, axes = plt.subplots(1, 2, figsize=(6.3, 3.2))
# 明确类型名称及顺序，使转移矩阵、初始向量、输出各列和图例保持对应。
names = ["农地", "建设用地", "其他用地"]
# errorbar 的 yerr 接收相对中心点的下/上偏移量，而不是直接接收上下端点。
# 误差棒以中位数为中心，允许上下不对称，表示边际情景分位区间而非标准误或置信区间。
axes[0].errorbar(np.arange(3), mid, yerr=[mid - lo, hi - mid],
                fmt="o", capsize=5)
# 叉号为固定 P 的十步期望，占比经非线性的 Q¹⁰ 变换后，情景中位数不必与叉号重合。
axes[0].scatter(np.arange(3), deterministic, marker="x", label="固定矩阵 P")
# 横轴只表示三类土地标签；纵轴是第十步期望占比，不是面积计数或抽样频率。
axes[0].set(title="第十步的边际情景区间",
            xticks=range(3), xticklabels=names,
            ylabel="期望面积占比", ylim=(0, 0.7))
# 图例标明固定矩阵点值，其余圆点与误差棒来自情景的中位数和分位数。
axes[0].legend(fontsize=8)
# 依次读取 pred 的三列，各列都是跨 2000 份情景的一维边际结果。
for j, name in enumerate(names):
    # 叠加三个半透明直方图；纵轴为情景数，未做概率密度归一化。
    axes[1].hist(pred[:, j], bins=30, alpha=0.4, label=name)
# 横轴强调是条件期望占比，右图宽度由人为设定的 狄利克雷 情景扰动决定。
axes[1].set(title="情景中的期望占比",
            xlabel="第十步期望面积占比", ylabel="情景频数")
# 添加类型图例并调整整体边距，避免长类型名称挤压坐标标签。
axes[1].legend(fontsize=8)
fig.tight_layout(pad=0.8, w_pad=1.2)
# 以 180 dpi 保存图并关闭对象，结束这个不依赖外部数据的应用示例。
# bbox_inches="tight" 防止面板标题和类别名称被裁切；教材直接嵌入此 PNG。
fig.savefig("application_land.png", dpi=180, bbox_inches="tight")
plt.close(fig)