注：书中对代码的讲解并不详细，本文对很多细节做了详细注释。另外，书上的源代码是在Jupyter Notebook上运行的，较为分散，本文将代码集中起来，并加以完善，全部用vscode在python 3.9.18下测试通过，同时对于书上部分章节也做了整合。

Chapter8 Recurrent Neural Networks

8.1 Sequence Models

用$x_t$表示在时间步（time step）$t\in {\mathbb{Z}}^{+}$时观察到的股票价格$x_t$。注意，$t$对于本文中的序列通常是离散的，并在整数或其子集上变化。假设一个交易员通过以下途径预测$x_t$：
$$x_t\sim P(x_t\mid x_{t-1},\dots ,x_1).$$

8.1.1 Autoregressive Models

为了实现这个预测，交易员可以使用回归模型，但输入数量以及$x_{t-1},\dots ,x_1$本身因$t$而异，有效估计$P(x_t\mid x_{t-1},\dots ,x_1)$展开归结为以下两种策略。
第一种策略是，假设在现实情况下相当长的序列$x_{t-1},\dots ,x_1$可能是不必要的，因此我们只需要满足某个长度为$\tau$的时间跨度，即使用观测序列$x_{t-1},\dots ,x_{t-\tau }$，好处就是参数的数量总是不变的（至少在$t>\tau$时如此），这种模型被称为自回归模型（autoregressive models），因为它们是对自己执行回归。

第二种策略如下图所示，指保留一些对过去观测的总结$h_t$，并同时更新预测${\hat{x}}_t$和总结$h_t$。这就产生了基于${\hat{x}}_t=P(x_t\mid h_t)$估计$x_t$，以及公式$h_t=g(h_{t-1},x_{t-1})$更新的模型。由于$h_t$从未被观测到，这类模型也被称为隐变量自回归模型（latent autoregressive models）。

一个常见的假设是虽然特定值$x_t$可能会改变，但是序列本身的动力不会改变。这样的假设是合理的，因为新的动力一定受新的数据影响，而我们不可能用目前所掌握的数据来预测新的动力。称不变的动力学静止的（stationary），因此整个序列的估计值都将通过以下的方式获得：

$$P(x_1,\dots ,x_T)=\prod _{t=1}^{T}P(x_t\mid x_{t-1},\dots ,x_1).$$

8.1.2 Markov Models

当$P(x_t\mid x_{t-1},\dots ,x_1)=P(x_t\mid x_{t-\tau },\dots ,x_{t-1})$时，称序列满足马尔可夫条件（Markov condition）。特别地，如果$\tau =1$，得到一阶马尔可夫模型（first-order Markov model），
$$P(x_1,\dots ,x_T)=\prod _{t=1}^{T}P(x_t\mid x_{t-1})\text{ , }P(x_1\mid x_0)=P(x_1).$$

当假设$x_t$仅是离散值时，使用动态规划可以沿着马尔可夫链精确地计算结果。例如，利用$P(x_{t+1}\mid x_t,x_{t-1})=P(x_{t+1}\mid x_t)$可以高效地计算$P(x_{t+1}\mid x_{t-1})$：

$$\begin{array}{rl}P(x_{t+1}\mid x_{t-1})& =\frac{{\sum }_{x_t}P(x_{t+1},x_t,x_{t-1})}{P(x_{t-1})}\\ & =\frac{{\sum }_{x_t}P(x_{t+1}\mid x_t,x_{t-1})P(x_t,x_{t-1})}{P(x_{t-1})}\\ & =\sum _{x_t}P(x_{t+1}\mid x_t)P(x_t\mid x_{t-1})\end{array}$$

原则上，将$P(x_1,\dots ,x_T)$倒序展开也没什么问题，因为基于条件概率公式有：

$$P(x_1,\dots ,x_T)=\prod _{t=T}^{1}P(x_t\mid x_{t+1},\dots ,x_T).$$
然而，在许多情况下，数据在时间上是前进的，显然未来的事件不能影响过去。在某些情况下，对于某些可加性噪声$ϵ$，我们可以找到$x_{t+1}=f(x_t)+ϵ$，反之则不行。

import torch
from torch import nn
from d2l import torch as d2l
import matplotlib.pyplot as pltT = 1000  # 总共产生1000个点
time = torch.arange(1, T + 1, dtype=torch.float32)
x = torch.sin(0.01 * time) + torch.normal(0, 0.2, (T,))
d2l.plot(time, [x], 'time', 'x', xlim=[1, 1000], figsize=(6, 3))
plt.show()

#将这个序列转换为模型的特征－标签（feature-label）对
tau = 4
features = torch.zeros((T - tau, tau))#这比提供的样本少了tau个，因为没有足够的历史记录来描述前tau个样本。
for i in range(tau):features[:, i] = x[i: T - tau + i]
labels = x[tau:].reshape((-1, 1))batch_size, n_train = 16, 600
# 只有前n_train个样本用于训练
train_iter = d2l.load_array((features[:n_train], labels[:n_train]),batch_size, is_train=True)# 初始化网络权重的函数
def init_weights(m):if type(m) == nn.Linear:nn.init.xavier_uniform_(m.weight)# 一个简单的多层感知机
def get_net():net = nn.Sequential(nn.Linear(4, 10),nn.ReLU(),nn.Linear(10, 1))net.apply(init_weights)return net# 平方损失。注意：MSELoss计算平方误差时不带系数1/2
loss = nn.MSELoss(reduction='none')def train(net, train_iter, loss, epochs, lr):trainer = torch.optim.Adam(net.parameters(), lr)for epoch in range(epochs):for X, y in train_iter:trainer.zero_grad()l = loss(net(X), y)l.sum().backward()trainer.step()print(f'epoch {epoch + 1}, 'f'loss: {d2l.evaluate_loss(net, train_iter, loss):f}')net = get_net()
train(net, train_iter, loss, 5, 0.01)#单步预测
onestep_preds = net(features)
d2l.plot([time, time[tau:]],[x.detach().numpy(), onestep_preds.detach().numpy()], 'time','x', legend=['data', '1-step preds'], xlim=[1, 1000],figsize=(6, 3))
plt.show()

即使预测的时间步超过了$600+4$（n_train + tau），单步预测效果不错，但如果数据观察序列的时间步只到$604$，我们需要一步一步地向前迈进：
$$\begin{gathered} {\hat{x}}_{605}=f(x_{601},x_{602},x_{603},x_{604}), \\ {\hat{x}}_{606}=f(x_{602},x_{603},x_{604},{\hat{x}}_{605}), \\ {\hat{x}}_{607}=f(x_{603},x_{604},{\hat{x}}_{605},{\hat{x}}_{606}), \\ {\hat{x}}_{608}=f(x_{604},{\hat{x}}_{605},{\hat{x}}_{606},{\hat{x}}_{607}), \\ {\hat{x}}_{609}=f({\hat{x}}_{605},{\hat{x}}_{606},{\hat{x}}_{607},{\hat{x}}_{608}), \\ \dots \end{gathered}$$
通常，对于直到$x_t$的观测序列，其在时间步$t+k$处的预测输出${\hat{x}}_{t+k}$称为$k$步预测（$k$-step-ahead-prediction）。

#多步预测
multistep_preds = torch.zeros(T)
multistep_preds[: n_train + tau] = x[: n_train + tau]
for i in range(n_train + tau, T):multistep_preds[i] = net(multistep_preds[i - tau:i].reshape((1, -1)))d2l.plot([time, time[tau:], time[n_train + tau:]],[x.detach().numpy(), onestep_preds.detach().numpy(),multistep_preds[n_train + tau:].detach().numpy()], 'time','x', legend=['data', '1-step preds', 'multistep preds'],xlim=[1, 1000], figsize=(6, 3))
plt.show()

如上面例子所示，绿线的预测经过几个步骤之后就会衰减到一个常数，原因在于错误的累积：假设在步骤$1$之后，我们积累了一些错误${ϵ}_1=\overline{ϵ}$。于是，步骤$2$的输入被扰动了${ϵ}_1$，结果积累的误差是${ϵ}_2=\overline{ϵ}+c{ϵ}_1$（$c$为常数），后面依此类推。基于$k=1,4,16,64$，通过对整个序列预测的计算，让我们更仔细地看一下$k$步预测的困难。

max_steps = 64
features = torch.zeros((T - tau - max_steps + 1, tau + max_steps))
# 列i（i<tau）是来自x的观测，其时间步从（i）到（i+T-tau-max_steps+1）
for i in range(tau):features[:, i] = x[i: i + T - tau - max_steps + 1]# 列i（i>=tau）是来自（i-tau+1）步的预测，其时间步从（i）到（i+T-tau-max_steps+1）
for i in range(tau, tau + max_steps):features[:, i] = net(features[:, i - tau:i]).reshape(-1)steps = (1, 4, 16, 64)
d2l.plot([time[tau + i - 1: T - max_steps + i] for i in steps],[features[:, tau + i - 1].detach().numpy() for i in steps], 'time', 'x',legend=[f'{i}-step preds' for i in steps], xlim=[5, 1000],figsize=(6, 3))
plt.show()