用Python实现时间序列模型实战——Day 24: 时间序列中的贝叶斯方法

一、学习内容

1. 贝叶斯时间序列分析的基本概念

贝叶斯方法基于贝叶斯统计，通过对数据的先验分布和似然函数进行推断，更新为后验分布。贝叶斯时间序列分析使用贝叶斯推断处理时间序列中的不确定性，适合处理复杂、不确定性高的时间序列问题。

贝叶斯时间序列模型能够自然地融合不同的信息源，并对不确定性进行定量描述。贝叶斯结构时间序列模型 (BSTS) 是一种常见的贝叶斯时间序列模型，它结合了状态空间模型和贝叶斯推断来进行建模和预测。

2. 基于贝叶斯推断的时间序列模型

贝叶斯时间序列模型的基本过程是通过先验分布、似然函数和后验分布来进行推断。常见的贝叶斯时间序列模型包括：

BSTS（贝叶斯结构时间序列）：适合处理带有趋势、季节性、回归等特征的时间序列。
贝叶斯动态线性模型（Bayesian Dynamic Linear Model, DLM）：适用于随时间变化的线性动态系统。

3. 时间序列中的不确定性处理与预测

贝叶斯方法的一个优势是能够很好地处理不确定性，并通过后验分布的预测结果，给出预测的置信区间。在时间序列预测中，可以通过预测的后验分布来量化预测的不确定性。

二、实战案例

面是使用 PyTorch 和 Pyro 构建贝叶斯时间序列模型的代码。

import torch
import torch.nn as nn
import pyro
import pyro.distributions as dist
from pyro.infer import SVI, Trace_ELBO
from pyro.optim import Adam
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.preprocessing import StandardScaler# 1. 数据加载与预处理
url = 'https://raw.githubusercontent.com/jbrownlee/Datasets/master/airline-passengers.csv'
data = pd.read_csv(url, header=0, parse_dates=['Month'], index_col='Month')# 数据标准化处理
scaler = StandardScaler()
scaled_data = scaler.fit_transform(data[['Passengers']])# 生成滞后特征
def create_lagged_features(data, lag=1):X, y = [], []for i in range(len(data) - lag):X.append(data[i:i + lag])y.append(data[i + lag])return np.array(X), np.array(y)# 使用滞后步长为12个月的数据
lag = 12
X, y = create_lagged_features(scaled_data, lag=lag)# 划分训练集和测试集
train_size = int(len(X) * 0.8)
X_train, X_test = torch.tensor(X[:train_size], dtype=torch.float32), torch.tensor(X[train_size:], dtype=torch.float32)
y_train, y_test = torch.tensor(y[:train_size], dtype=torch.float32), torch.tensor(y[train_size:], dtype=torch.float32)# 2. PyTorch 模型
class BayesianLSTM(nn.Module):def __init__(self, input_size, hidden_size):super(BayesianLSTM, self).__init__()self.lstm = nn.LSTM(input_size, hidden_size, batch_first=True)self.fc = nn.Linear(hidden_size, 1)def forward(self, x):lstm_out, _ = self.lstm(x)output = self.fc(lstm_out[:, -1, :])return output# 3. Pyro 模型定义
def model(x_data, y_data=None):# 定义先验分布hidden_size = 50input_size = x_data.shape[-1]# 权重的先验分布lstm_weight_prior = dist.Normal(torch.zeros(hidden_size, input_size), torch.ones(hidden_size, input_size))lstm_bias_prior = dist.Normal(torch.zeros(hidden_size), torch.ones(hidden_size))# 全连接层的权重与偏置的先验fc_weight_prior = dist.Normal(torch.zeros(1, hidden_size), torch.ones(1, hidden_size))fc_bias_prior = dist.Normal(torch.zeros(1), torch.ones(1))# 定义 Pyro 参数priors = {'lstm.weight_ih_l0': lstm_weight_prior,'lstm.bias_ih_l0': lstm_bias_prior,'fc.weight': fc_weight_prior,'fc.bias': fc_bias_prior}# 使用 Pyro 模型bayesian_lstm = pyro.module('bayesian_lstm', BayesianLSTM(input_size, hidden_size))# 预测prediction = bayesian_lstm(x_data)# 似然函数with pyro.plate('data', len(x_data)):likelihood = dist.Normal(prediction, torch.tensor(0.1))obs = likelihood.log_prob(y_data) if y_data is not None else likelihoodreturn obs# 4. 损失函数和优化器
def guide(x_data, y_data=None):# 引导函数用于推断模型参数pyro.module('bayesian_lstm', BayesianLSTM(input_size=X_train.shape[-1], hidden_size=50))# 使用 SVI 进行推断
optim = Adam({"lr": 0.01})
svi = SVI(model, guide, optim, loss=Trace_ELBO())# 5. 模型训练
num_iterations = 1000
losses = []
for epoch in range(num_iterations):loss = svi.step(X_train, y_train)losses.append(loss)if epoch % 100 == 0:print(f"Epoch {epoch}, Loss: {loss}")# 6. 模型预测
bayesian_lstm = BayesianLSTM(input_size=X_train.shape[-1], hidden_size=50)
bayesian_lstm.eval()
y_pred = bayesian_lstm(X_test).detach().numpy()# 还原到原始数据尺度
y_test_rescaled = scaler.inverse_transform(y_test)
y_pred_rescaled = scaler.inverse_transform(y_pred)# 7. 结果可视化
plt.figure(figsize=(12, 6))
plt.plot(data.index[-len(y_test):], y_test_rescaled, label='Actual Passengers')
plt.plot(data.index[-len(y_test):], y_pred_rescaled, label='Predicted Passengers', color='green')
plt.title('Bayesian Time Series Predictions with Pyro')
plt.xlabel('Date')
plt.ylabel('Number of Passengers')
plt.legend()
plt.grid(True)
plt.show()