DeepAR自回归概率预测采样

可视化演示

\u901a\u8fc7\u4e09\u9636\u6bb5\u4ea4\u4e92\u5c55\u793aDeepAR\u81ea\u56de\u5f52\u7f51\u7edc\u67b6\u6784\u3001\u8499\u7279\u5361\u6d1b\u91c7\u6837\u591a\u8def\u5f84\u4e0eCopula\u8054\u5408\u5206\u5e03\u5efa\u6a21\u3002

自回归原理

小林在工作中注意到一个朴素的规律：天气是有“记忆”的。昨天26度，今天大概率在24到28度之间；昨天下了暴雨，今天地面积水还没干，空气湿度依然很大。换句话说，未来的天气高度依赖于最近的过去。这种“用自己过去的值来预测自己下一步”的思路，在统计学里被称为“自回归”（Auto-Regressive）——DeepAR模型的名字正是由此而来。

DeepAR（Deep Autoregressive）是亚马逊在2017年提出的概率预测模型。它的核心构件是一个LSTM循环神经网络：先把最近$w$个时间步的观测值$z_{t-w+1}, \ldots, z_t$作为“历史窗口”送入网络，LSTM将它们逐时刻编码成一个隐状态$h_t$——你可以把这个隐状态理解成小林大脑里对“这几天天气走势”形成的综合直觉。有了$h_t$，模型就能输出下一个时刻的“条件概率分布”：在已知过去的前提下，明天最可能是什么温度、波动有多大。

$$\displaystyle p(z_{t+1} \mid z_{t}, z_{t-1}, \ldots, z_{t-w+1}) = \mathcal{N}\left(z_{t+1}; \mu_{t+1}(h_t), \sigma_{t+1}^2(h_t)\right)$$

1. 条件概率符号：竖线左边$z_{t+1}$是我们要预测的明天温度，竖线右边是已知的历史观测值，读作“在已知历史的条件下，明天的概率分布”。

$$p(z_{t+1} \mid z_{t}, z_{t-1}, \ldots, z_{t-w+1})$$

2. 历史窗口：最近$w$个时间步的观测值，相当于小林翻阅的“过去一周温度记录”——窗口越长，模型掌握的背景信息越丰富。

$$z_{t}, z_{t-1}, \ldots, z_{t-w+1}$$

3. 高斯分布符号：钟形概率密度函数，由均值和方差两个参数完全决定，是概率预测最常用的分布假设。

$$\mathcal{N}\left(z_{t+1}; \mu, \sigma^2\right)$$

4. 均值函数：网络从隐状态$h_t$中解码出的“最可能温度”，是小林预报里报给市民的那个数字。

$$\mu_{t+1}(h_t)$$

5. 方差函数：网络解码出的“不确定程度”，方差越大，曲线越矮胖，说明小林的底气越不足。

$$\sigma_{t+1}^2(h_t)$$

有了条件分布的表达方式，小林就明白了一个道理：所谓“自回归”，不是简单地复制历史，而是从历史的脉络中提取规律，再把它外推到未来。这就像他看过去七天的等温线图，心里默默推演“明天这股冷空气大概走到哪儿了”——每一步预测都踩在历史的基础上，滚动向前。

概率参数化

第二天，小林尝试向市民发布“明天26度”的预报。老周看完直摇头：“26度？你知道明天真正是26度的概率有多大吗？不到10%。你报的'26度'几乎必然出错。”小林愣住了。老周接着说：“你得这样报：明天最可能是26度，但大概率落在22到30度之间。你得把'大概'两个字报出来。”

这句话点醒了小林：预报的核心不是猜一个最可能的数字，而是刻画一整条“概率曲线”的形状。DeepAR正是这么做的——它不在最后一层输出一个预测值，而是输出两个值：均值$\mu$和方差$\sigma^2$。均值告诉人们“最可能在哪”，方差告诉人们“有多不确定”。为了让方差始终为正数，网络实际输出的是“对数方差”$\log\sigma^2$，再通过指数运算还原。而训练时使用的目标函数，是负对数似然（Negative Log-Likelihood，简称NLL）损失——本质上是在说：请给我一个分布，让“真实发生的事件”在这个分布下的概率尽量高。

$$\displaystyle \mathcal{L}_{\mathrm{NLL}} = \frac{1}{T}\sum_{t=1}^{T}\left[\frac{(z_{t+1}-\mu_{t+1})^2}{2\sigma_{t+1}^2} + \frac{1}{2}\log(2\pi\sigma_{t+1}^2)\right]$$

1. 求和与平均：把所有时间步的损失加总后取平均，这样不同长度序列的损失可以公平比较。

$$\frac{1}{T}\sum_{t=1}^{T}$$

2. 均方误差项：真实值偏离预测均值的平方，再除以$2\sigma^2$。方差越小、偏差越大，惩罚就越重——就像小林把区间报得很窄，结果第二天温度跌出区间，就要挨批评。

$$\frac{(z_{t+1}-\mu_{t+1})^2}{2\sigma_{t+1}^2}$$

3. 方差惩罚项：$\log(2\pi\sigma^2)$随方差增大而增大，防止模型偷懒把方差调到无穷大来逃避第一项惩罚。它提醒小林：区间不能无限宽，否则预报就失去了意义。

$$\frac{1}{2}\log(2\pi\sigma_{t+1}^2)$$

4. 对数方差的由来：因为$\log\sigma^2$可以是任意实数，网络输出它之后再指数化，就保证了方差永远为正——这是数值稳定性的第一道保险。

$$\sigma^2 = \exp(\log\sigma^2) > 0$$

小林恍然大悟：好的预报员不会为了“准”而把话说死，而是诚实地给出“均值加波动”的完整描述。DeepAR和优秀预报员的共同点是：它们优化的不是“预测值是否等于真实值”，而是“预测的整个分布是否与真实世界相符”。

蒙特卡洛采样

模型训练好了，小林面临一个新问题：DeepAR给出的只是“明天温度服从某个分布”——但市长要的是未来30天的温度曲线，不是一个分布表达式。小林问老周：“我总不能把数学公式贴在官网上吧？”老周笑了：“你不会电脑模拟吗？让电脑按你说的分布，随机生成一万条'可能未来'，把这一万条画出来，就是最好的预报图。”

这就是蒙特卡洛采样（Monte Carlo Sampling）的思想。DeepAR在预测阶段，从最后一步的观测值开始，重复以下操作：先从当前步预测的条件分布中随机抽一个值，把这个“虚构的观测值”当作下一步的输入，再预测下一步的分布，再抽样……如此反复，就得到一条完整的未来路径。重复很多次，就得到一束从“现在”出发、逐渐散开的路径。路径越分散，说明不确定性越大——这正是“预测越远越不可靠”的直观体现。

$$\displaystyle \tilde{z}_{t+1}^{(i)} = \mu_{t+1}^{(i)} + \sigma_{t+1}^{(i)} \cdot \varepsilon^{(i)},\quad \varepsilon^{(i)} \sim \mathcal{N}(0,1)$$

1. 采样路径值：上标$(i)$表示第$i$条采样路径在$t+1$时刻的值。一千次采样，就是一千条“平行宇宙”里的明天。

$$\tilde{z}_{t+1}^{(i)}$$

2. 均值项：该路径在当前时刻的分布中心，来自网络输出的$\mu_{t+1}^{(i)}$，是所有采样路径共享的“主心骨”。

$$\mu_{t+1}^{(i)}$$

3. 随机扰动项：标准差$\sigma_{t+1}^{(i)}$乘以标准正态噪声$\varepsilon^{(i)}$，给路径加入合理的随机波动——波动幅度恰好等于模型估计的不确定性。

$$\sigma_{t+1}^{(i)} \cdot \varepsilon^{(i)}$$

4. 标准正态噪声：$\varepsilon^{(i)}$从均值为0、方差为1的正态分布中抽取。将这个噪声乘上$\sigma$、加上$\mu$，就得到了服从$\mathcal{N}(\mu,\sigma^2)$的样本。

$$\varepsilon^{(i)} \sim \mathcal{N}(0,1)$$

当小林把一万条模拟路径画在同一张图上时，他被深深震撼了：前3天路径还紧紧抱成一束，像秋天麦田里整齐的稻穗；到了第10天，路径散成一把巨大的扇子。他对老周说：“我终于明白什么叫'不确定性随时间传播'了——不是我不想报准，是这个世界本身在变得越来越不确定。”

Copula建模

冷锋过境的夜晚，小林发现了一个麻烦：他手上有东城、西城、郊区三份独立的预报，但三份预报拼在一起却“打架”——东城说降温8度，西城说只降2度，郊区甚至说升温。老周连夜打来电话：“三个区离得这么近，同一个冷锋扫过来，怎么可能一个暴跌一个暴涨？你要报的是'联合天气'，不是三张各说各话的纸。”

这就引出了概率预测中的一个经典难题：如何让多个相关的序列共享同一个“天气剧本”。Copula（连接函数）就是答案。简单说，Copula是一管“相关性胶水”：它把每个变量自己的“边缘分布”（比如东城温度服从什么分布、西城温度服从什么分布）和变量之间的“依赖结构”（比如东城降温时西城同时降温的概率有多大）分开建模。Sklar定理告诉我们：任意一个联合分布，都可以拆成一个Copula函数加上各变量的边缘分布。

$$\displaystyle C_{\theta}(u_1,\ldots,u_d) = \exp\left(-\left(\sum_{j=1}^{d}(-\log u_j)^{\theta}\right)^{1/\theta}\right),\quad \theta \geq 1$$

1. 概率积分变换：$u_j = F_j(z_j)$把第$j$个变量的观测值映射到$(0,1)$区间，相当于问“这个值在自身分布里排第几”。

$$u_j = F_j(z_j)$$

2. 极端程度求和：$-\log u_j$把接近1的$u_j$放大成很大的正数，越是极端的事件贡献越大；对所有变量求和后，$1/\theta$次方将总和折算回合适的尺度。

$$\left(\sum_{j=1}^{d}(-\log u_j)^{\theta}\right)^{1/\theta}$$

3. 相关强度参数：$\theta$控制着变量间的相关强度。$\theta=1$时各变量相互独立，$\theta$越大，变量越倾向于同涨同跌。

$$\theta \geq 1$$

4. 指数包装：最外层的$\exp$保证Copula的函数值落在合法的概率区间内，相当于给整管“胶水”套上了一个合规的容器。

$$\exp(-(\cdot))$$

小林把三份独立预报的边缘分布输入Copula，重新采样后，三份预报终于“步调一致”了：冷锋来时，三区的预测曲线一起下探，只是幅度略有差异。他感慨地说：“Copula不是把结果拼在一起，而是让它们共享同一个剧本——这样市民看到的才是真实世界的联动，而不是三个互不相干的故事。”

损失评估

一个月后，小林信心满满地发布了一份“未来7天：90%概率温度落在15到28度之间”的预报。结果一周后复盘，老周把实际数据拍在桌上：“7天里有5天落在区间里，覆盖率只有71%。你嘴上说着90%的把握，实际却没那个底气——这叫'校准失败'。你的区间报宽一点，覆盖率自然上去；但报得无限宽，预报就变成了废话。关键是怎么平衡。”

概率预测的评估，从来不是单一指标说了算。负对数似然（NLL）可以比较两个模型谁的分布更贴近真实世界，但它不够直观。更常用的评估工具是分位数损失（pinball loss）和区间覆盖率（PICP，即真实值落在预测区间内的比例）。分位数损失专为分位数回归设计：给定目标分位数$\tau$（比如0.9），模型输出的预测值应该是$\tau$分位数——这意味着真实值以$\tau$的概率低于它、以$1-\tau$的概率高于它。这个损失不对称地惩罚“预测过高”和“预测过低”。

$$\displaystyle \rho_{\tau}(u) = u \cdot \left(\tau - \mathbb{1}_{u<0}\right),\quad u = y - \hat{q}_{\tau}$$

1. 目标分位数：$\tau$是我们想要的分位水平，如$\tau=0.9$表示“预测上限要盖住90%的情况”。

$$\tau \in (0,1)$$

2. 预测分位数：$\hat{q}_{\tau}$是模型输出的第$\tau$分位数估计，也就是预测区间的边界。

$$\hat{q}_{\tau}$$

3. 预测残差：$u$是真实值超出预测分位的部分，正数表示真实值更高（模型低估），负数表示真实值更低（模型高估）。

$$u = y - \hat{q}_{\tau}$$

4. 指示函数与不对称权重：$\mathbb{1}_{u<0}$在$u<0$时取1，此时惩罚系数为$\tau-1$（负权重乘以负数，结果为正）；$u \geq 0$时惩罚系数为$\tau$。这就实现了“一边罚得重、一边罚得轻”的不对称设计。

$$\mathbb{1}_{u<0} = \begin{cases} 1, & u<0 \\ 0, & u \geq 0 \end{cases}$$

老周看着小林画出的分位损失曲线说：“记住，预测区间不能只看宽度，也不能只看覆盖率。宽度太窄、覆盖率不够，说明你盲目自信；宽度太宽、覆盖率超标，说明你在'用废话避险'。只有用分位数损失反复校验，才能让预报既锋利又诚实。”小林点点头——概率预测的评估，本质上是在“自信”与“准确”之间走钢丝。

工程实现

系统上线第三天，凌晨三点，监控警报把小林从睡梦中吵醒。屏幕上是一串刺眼的红色NaN——某个传感器回传了一个异常值，经过LSTM的层层传递，梯度爆炸，整个预报系统崩溃了。小林手忙脚乱地重启服务，晨间预报差点开天窗。第二天他红着眼问老周：“模型在论文里好好的，怎么一上生产就崩溃？”

老周摆摆手：“做研究和做工程的差别就在这里。概率预测模型的工程实现，一半的功夫都在'防坍塌'上。”他带小林逐项检查：第一，方差必须恒为正——网络输出原始值$o$后，用softplus函数$\log(1+e^o)$映射，任何实数输入都能变成正数，而且不会像指数函数那样容易溢出；第二，LSTM的梯度容易爆炸，训练时要加梯度裁剪；第三，推理阶段要一次性批量采样几千条路径，用矩阵运算代替逐条循环，才能满足每天清晨的出报时效；第四，窗口长度、隐藏层维度、学习率这些超参数，需要用网格搜索或贝叶斯优化仔细调。

$$\displaystyle \sigma_{t+1} = \log\left(1 + \exp(o_{t+1})\right)$$

1. 网络原始输出：$o_{t+1}$是LSTM解码器的原始输出，它可以取任意实数——正是这种“任意性”可能导致标准差为负，产生数学上的荒谬。

$$o_{t+1} \in \mathbb{R}$$

2. 指数映射：$\exp(o)$把实数映射到正数域，保证基本非负，但纯指数在$o$较大时会急剧膨胀，所以需要再加一层保护。

$$\exp(o_{t+1})$$

3. 对数包装：$\log(1+\cdot)$是softplus的精髓：当$o$为很大的负数时，结果趋近于0但绝不等于0；当$o$为很大的正数时，结果近似于$o$本身。曲线平滑可导，梯度稳定。

$$\log(1+\exp(o_{t+1}))$$

小林连夜给系统加上了三道“护栏”：softplus保方差为正、梯度裁剪防爆炸、批量采样提速度。第二天清晨，预报系统准时把全市的联合概率预报推送到了市民手机上。他看着屏幕上平稳运行的健康指标松了口气——工程化没有那么多浪漫，它的本质是预判每一种崩溃的可能，并在崩溃发生之前，用一行行朴素的代码把它拦住。

AI 应用场景

Python 演示

import numpy as np # 生成带噪声的正弦波序列np.random.seed(42)T = 300t = np.arange(T) / 10.0z = np.sin(t) * 2.0 + np.random.randn(T) * 0.2 # 构造滑动窗口样本:用过去 WINDOW 个点预测下一点WINDOW = 8X, y = [], []for i in range(WINDOW, T - 1):    X.append(z[i - WINDOW:i])    y.append(z[i])X = np.array(X).reshape(-1, WINDOW)y = np.array(y) # 两层网络:输出均值 mu 与对数方差 log_varrng = np.random.default_rng(0)W1 = rng.standard_normal((WINDOW, 32)) * 0.5b1 = np.zeros(32)W2_mu = rng.standard_normal((32, 1)) * 0.5W2_lv = rng.standard_normal((32, 1)) * 0.1lr = 0.01 def forward(x):    h = np.tanh(x @ W1 + b1)    mu = h @ W2_mu    log_var = h @ W2_lv    return mu.flatten(), log_var.flatten() def nll_loss(y, mu, log_var):    var = np.exp(log_var)    return np.mean((y - mu) ** 2 / (2 * var) + 0.5 * log_var) # 手动反向传播训练 200 轮for _ in range(200):    mu, lv = forward(X)    var = np.exp(lv)    dmu = -(y - mu) / var    dlv = (y - mu) ** 2 / (2 * var) - 0.5    h = np.tanh(X @ W1 + b1)    W2_mu -= lr * (h.T @ dmu.reshape(-1, 1)) / len(X)    W2_lv -= lr * (h.T @ dlv.reshape(-1, 1)) / len(X)    db = (dmu.reshape(-1, 1) @ W2_mu.T + dlv.reshape(-1, 1) @ W2_lv.T) * (1 - h ** 2)    W1 -= lr * (X.T @ db) / len(X)    b1 -= lr * db.mean(axis=0) * 0.1 mu, lv = forward(X)print('训练后 NLL:', nll_loss(y, mu, lv)) # 蒙特卡洛多步采样:从历史末端滚动生成未来路径def mc_forecast(history, steps=10, n_samples=100):    samples = []    for _ in range(n_samples):        path = list(history)        for _ in range(steps):            x = np.array(path[-WINDOW:]).reshape(1, -1)            m, lv = forward(x)            path.append(m[0] + np.exp(lv[0] * 0.5) * np.random.randn())        samples.append(path[-steps:])    return np.array(samples) preds = mc_forecast(z[:80], steps=10, n_samples=200)lo = np.percentile(preds, 10, axis=0)hi = np.percentile(preds, 90, axis=0)print('90% 区间示例:')for i in range(3):    print('  步', i + 1, ':', round(lo[i], 2), '~', round(hi[i], 2))

1. z = np.sin(t) * 2.0 + np.random.randn(T) * 0.2生成带噪声的周期序列，模拟“有规律但不确定”的时间信号
2. X.append(z[i - WINDOW:i])构造自回归窗口：用过去8个点预测下一个点
3. W2_mu = rng.standard_normal((32, 1)) * 0.5均值输出头的权重，负责产生“最可能的预测”
4. W2_lv = rng.standard_normal((32, 1)) * 0.1对数方差输出头的权重，初始化更小，避免起始方差过大
5. var = np.exp(lv)将对数方差还原为方差，保证方差恒为正
6. (y - mu) ** 2 / (2 * var) + 0.5 * log_var负对数似然的两项：均方误差按方差缩放 + 方差自身惩罚
7. dlv = (y - mu) ** 2 / (2 * var) - 0.5对对数方差的梯度，让模型学会“何时该自信、何时该保守”
8. m[0] + np.exp(lv[0] * 0.5) * np.random.randn()重参数化采样：均值加标准差乘标准正态噪声
9. np.percentile(preds, 10, axis=0)对200条路径逐时间步取10%分位，得到区间下界

常见误区

面试题