时间序列预测方法¶

import numpy as np
import pandas as pd
import torch
import torch.nn as nn

# ===== 1. ARIMA 统计方法（用 statsmodels）=====
# pip install statsmodels
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.stattools import adfuller

# 模拟微生物丰度时间序列（周采样，52周）
np.random.seed(42)
microbiome_abundance = np.random.randn(52).cumsum() + 10  # 随机游走模拟

# ADF 检验平稳性
result = adfuller(microbiome_abundance)
print(f'ADF 检验 p 值: {result[1]:.4f}')
# p < 0.05: 平稳；p > 0.05: 非平稳需差分

# 拟合 ARIMA 模型
# p=2: 自回归阶数；d=1: 差分阶数；q=1: 移动平均阶数
model = ARIMA(microbiome_abundance, order=(2, 1, 1))
fitted = model.fit()
print(fitted.summary())

# 预测未来 8 周
forecast = fitted.forecast(steps=8)
print(f'未来8周预测: {forecast.round(2)}')

# ===== 2. LSTM 深度学习方法 =====
def create_sequences(data, seq_len=10, pred_len=1):
    """
    把时间序列切成滑动窗口训练样本
    seq_len: 用多少历史步预测
    pred_len: 预测未来多少步
    """
    X, y = [], []
    for i in range(len(data) - seq_len - pred_len + 1):
        X.append(data[i:i+seq_len])          # 输入：过去 seq_len 步
        y.append(data[i+seq_len:i+seq_len+pred_len])  # 标签：未来 pred_len 步
    return np.array(X), np.array(y)

class LSTMForecaster(nn.Module):
    """LSTM 时间序列预测器"""
    def __init__(self, input_size=1, hidden_size=64, num_layers=2, pred_len=1):
        super().__init__()
        # LSTM 层：处理时间序列依赖
        self.lstm = nn.LSTM(
            input_size=input_size,    # 每个时间步的特征维度（单变量=1）
            hidden_size=hidden_size,  # 隐藏状态维度
            num_layers=num_layers,    # 堆叠 LSTM 层数
            batch_first=True,         # 输入格式: [batch, seq, features]
            dropout=0.2               # 层间 dropout
        )
        # 线性输出层：LSTM 输出 -> 预测值
        self.fc = nn.Linear(hidden_size, pred_len)

    def forward(self, x):
        # x 形状: [batch, seq_len, input_size]
        lstm_out, _ = self.lstm(x)      # lstm_out: [batch, seq_len, hidden_size]
        last_step = lstm_out[:, -1, :]  # 取最后一步的输出 [batch, hidden_size]
        prediction = self.fc(last_step) # 线性变换到预测长度 [batch, pred_len]
        return prediction

# 数据准备
seq_len = 10   # 用过去 10 周
pred_len = 1   # 预测下 1 周

X, y = create_sequences(microbiome_abundance, seq_len, pred_len)
X = torch.FloatTensor(X).unsqueeze(-1)  # [N, seq_len, 1]
y = torch.FloatTensor(y)                 # [N, pred_len]

# 数据归一化（时间序列必须！）
X_mean, X_std = X.mean(), X.std()
X = (X - X_mean) / X_std
y = (y - X_mean) / X_std

# 划分训练/测试集
train_size = int(0.8 * len(X))
X_train, X_test = X[:train_size], X[train_size:]
y_train, y_test = y[:train_size], y[train_size:]

# 训练
model = LSTMForecaster(input_size=1, hidden_size=32, num_layers=2, pred_len=pred_len)
optimizer = torch.optim.Adam(model.parameters(), lr=1e-3)
loss_fn = nn.MSELoss()

for epoch in range(100):
    model.train()
    pred = model(X_train)
    loss = loss_fn(pred, y_train)
    optimizer.zero_grad()
    loss.backward()
    optimizer.step()
    if epoch % 20 == 0:
        print(f'Epoch {epoch}, Loss: {loss.item():.4f}')

# 评估（反归一化）
model.eval()
with torch.no_grad():
    pred_test = model(X_test) * X_std + X_mean  # 反归一化
    true_test = y_test * X_std + X_mean
    mae = torch.abs(pred_test - true_test).mean()
    print(f'测试 MAE: {mae.item():.4f}')