Diebold-Mariano 检验 Python 实战：3步完成模型预测精度显著性对比

Diebold-Mariano 检验 Python 实战3步完成模型预测精度显著性对比在时间序列预测项目中我们常常需要比较不同模型的预测效果。比如ARIMA和LSTM模型在同一数据集上的表现孰优孰劣仅仅观察平均误差的差异可能不够严谨——这差异是真实的模型能力体现还是数据采样带来的偶然结果Diebold-Mariano检验简称DM检验正是为解决这一问题而生。本文将带你从零实现DM检验的Python代码并通过完整案例演示如何在实际项目中应用。无论你是数据分析师还是算法工程师都能快速掌握这套方法论为模型选择提供统计显著性依据。1. DM检验核心原理与Python函数封装DM检验的核心思想是通过统计检验判断两个模型的预测误差差异是否显著。其零假设是两个模型的预测精度相同若检验结果拒绝该假设则说明一个模型显著优于另一个。关键计算步骤计算损失差分序列对于两个模型的预测误差序列e1和e2长度均为T计算每个时间点的损失差分d [e1[i] - e2[i] for i in range(T)] # 假设使用MSE损失计算统计量DM统计量服从标准正态分布计算公式为import numpy as np d_mean np.mean(d) d_std np.std(d, ddof1) # 样本标准差 DM d_mean / (d_std / np.sqrt(T))结果解读查标准正态分布表得到p-valuefrom scipy import stats p_value 2 * (1 - stats.norm.cdf(abs(DM))) # 双侧检验完整Python函数实现def dm_test(e1, e2, loss_funcmse, h1, power2): Diebold-Mariano检验实现 参数: e1, e2: 两个模型的预测误差序列(长度相同) loss_func: 损失函数类型(mse, mae或自定义函数) h: 预测步长(用于调整自相关) power: 损失函数的幂次(当loss_funcmse时自动设为2) 返回: DM统计量, p-value assert len(e1) len(e2), 误差序列长度必须相同 T len(e1) # 计算损失差分 if loss_func mse: d [e1[i]**2 - e2[i]**2 for i in range(T)] elif loss_func mae: d [abs(e1[i]) - abs(e2[i]) for i in range(T)] else: d [loss_func(e1[i]) - loss_func(e2[i]) for i in range(T)] # 计算DM统计量(考虑自相关) d_mean np.mean(d) gamma [np.mean([(d[k] - d_mean) * (d[k - lag] - d_mean) for k in range(lag, T)]) for lag in range(1, h)] var_d np.var(d, ddof1) 2 * sum(gamma) DM d_mean / np.sqrt(var_d / T) # 计算p-value p_value 2 * (1 - stats.norm.cdf(abs(DM))) return DM, p_value提示实际应用中预测误差可能存在自相关特别是多步预测时。上述实现通过参数h考虑了这一点使用Newey-West方差估计器进行修正。2. 端到端应用案例ARIMA vs LSTM模型比较让我们通过一个完整案例演示DM检验的实际应用。假设我们有一个月度销售数据集分别用ARIMA和LSTM模型进行预测现在要比较两者的预测精度。数据准备与模型训练import pandas as pd from statsmodels.tsa.arima.model import ARIMA from tensorflow.keras.models import Sequential from tensorflow.keras.layers import LSTM, Dense # 加载数据 data pd.read_csv(sales.csv, parse_dates[date], index_coldate) train data.iloc[:-24] # 前N-24个月训练 test data.iloc[-24:] # 最后24个月测试 # ARIMA模型 arima ARIMA(train, order(1,1,1)).fit() arima_pred arima.forecast(24) # LSTM模型准备 def create_dataset(dataset, look_back1): X, Y [], [] for i in range(len(dataset)-look_back): X.append(dataset[i:(ilook_back)]) Y.append(dataset[ilook_back]) return np.array(X), np.array(Y) look_back 3 train_X, train_Y create_dataset(train.values, look_back) test_X, test_Y create_dataset(test.values, look_back) # LSTM模型训练 model Sequential() model.add(LSTM(50, input_shape(look_back, 1))) model.add(Dense(1)) model.compile(lossmean_squared_error, optimizeradam) model.fit(train_X, train_Y, epochs100, batch_size1, verbose0) lstm_pred model.predict(test_X)计算预测误差并执行DM检验# 计算预测误差 arima_errors test.values.flatten() - arima_pred lstm_errors test_Y.flatten() - lstm_pred.flatten() # 执行DM检验 dm_stat, p_value dm_test(arima_errors, lstm_errors, loss_funcmse) print(fDM统计量: {dm_stat:.4f}) print(fP-value: {p_value:.4f}) if p_value 0.05: if dm_stat 0: print(结论: ARIMA预测显著优于LSTM (p 0.05)) else: print(结论: LSTM预测显著优于ARIMA (p 0.05)) else: print(结论: 两个模型预测精度无显著差异)结果可视化分析import matplotlib.pyplot as plt plt.figure(figsize(12, 6)) plt.plot(test.index, test.values, label实际值, colorblack) plt.plot(test.index, arima_pred, labelARIMA预测, linestyle--) plt.plot(test.index[look_back:], lstm_pred, labelLSTM预测, linestyle:) plt.title(模型预测效果对比) plt.legend() plt.show() # 绘制误差分布 plt.figure(figsize(10, 5)) plt.boxplot([arima_errors, lstm_errors], labels[ARIMA误差, LSTM误差]) plt.title(预测误差分布对比) plt.show()3. 高级应用与注意事项3.1 不同损失函数的选择DM检验支持多种损失函数常见选择包括损失函数公式适用场景MSEe²强调大误差惩罚MAEeMAPEe/y# 使用MAE损失的DM检验 dm_mae, p_mae dm_test(arima_errors, lstm_errors, loss_funcmae) # 自定义损失函数示例 def huber_loss(e, delta1.0): return np.where(np.abs(e) delta, 0.5*e**2, delta*(np.abs(e)-0.5*delta)) dm_huber, p_huber dm_test(arima_errors, lstm_errors, loss_funchuber_loss)3.2 多步预测的调整当进行h步预测时误差通常呈现MA(h-1)自相关结构。这时需要设置合适的h参数使用Newey-West方差估计考虑小样本修正HLN检验# 12步超前预测的DM检验 dm_multi, p_multi dm_test(arima_errors, lstm_errors, h12) # HLN修正(小样本更准确) def hln_test(dm_stat, h, T): return dm_stat / np.sqrt((T 1 - 2*h h*(h-1)/T) / T) hln_stat hln_test(dm_multi, h12, Tlen(arima_errors)) hln_p 2 * (1 - stats.norm.cdf(abs(hln_stat)))3.3 常见问题解决方案问题1检验结果与误差均值矛盾可能原因误差分布不对称或存在极端值解决方案尝试不同损失函数或检查误差分布问题2p值边界显著如0.04-0.06建议增加测试集样本量或使用交叉验证多次检验问题3多重检验问题当比较多个模型时需要进行p值校正from statsmodels.stats.multitest import multipletests pvals [0.03, 0.05, 0.01] # 假设三个比较的p值 reject, adj_pvals, _, _ multipletests(pvals, methodholm)4. 工程实践中的优化技巧在实际项目中我们可以进一步优化DM检验的实现4.1 自动化检验流程def auto_dm_test(model_dict, test_data, loss_funcs[mse, mae]): 自动化多模型比较 model_dict: {模型名称: 预测值}的字典 test_data: 真实值 loss_funcs: 使用的损失函数列表 results [] models list(model_dict.keys()) for i in range(len(models)): for j in range(i1, len(models)): e1 test_data - model_dict[models[i]] e2 test_data - model_dict[models[j]] for loss in loss_funcs: dm, p dm_test(e1, e2, loss_funcloss) results.append({ model1: models[i], model2: models[j], loss_func: loss, DM_stat: dm, p_value: p, significance: p 0.05 }) return pd.DataFrame(results)4.2 结合交叉验证from sklearn.model_selection import TimeSeriesSplit def cv_dm_test(model1, model2, data, n_splits5): 交叉验证版DM检验 tscv TimeSeriesSplit(n_splitsn_splits) p_values [] for train_idx, test_idx in tscv.split(data): train, test data[train_idx], data[test_idx] # 训练模型并预测(此处简化实际需适配具体模型) pred1 model1.fit(train).predict(test) pred2 model2.fit(train).predict(test) _, p dm_test(test-pred1, test-pred2) p_values.append(p) # 组合各折p值(Fisher方法) combined_stat -2 * sum(np.log(p_values)) combined_p 1 - stats.chi2.cdf(combined_stat, df2*n_splits) return combined_p4.3 性能优化建议对于大规模时间序列可以使用numpy向量化计算对长序列采用滚动窗口检验并行化多模型比较# 向量化实现示例 def dm_test_vectorized(e1, e2, h1): d e1**2 - e2**2 # MSE d_mean np.mean(d) lags np.arange(1, h) gamma [np.sum((d[lag:] - d_mean) * (d[:-lag] - d_mean)) / len(d) for lag in lags] var_d np.var(d, ddof1) 2 * sum(gamma) DM d_mean / np.sqrt(var_d / len(d)) p 2 * (1 - stats.norm.cdf(abs(DM))) return DM, p在实际项目中我曾遇到一个有趣案例当比较Prophet和XGBoost在销售预测中的表现时DM检验结果随预测周期长度变化而反转。短期1-3天XGBoost显著更优但长期7天以上Prophet反而更好。这提醒我们模型比较需要放在具体应用场景下评估。