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

写代码的M教授

于 2024-09-18 08:00:00 发布

阅读量223

点赞数 6

分类专栏：时间序列模型文章标签： python 开发语言

本文链接：https://blog.csdn.net/qq_41698317/article/details/142304208

版权

时间序列模型专栏收录该内容

25 篇文章 1 订阅

订阅专栏

一、学习内容

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 = 50
    input_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 likelihood
        return 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()