sklearn线性与多项式回归如何计算预测置信区间

sklearn 没有内置直接输出预测置信/预测区间的功能，所有线性类模型（包括多项式回归，本质是特征转换后的线性回归）的80%预测区间都可以通过统计公式手动计算，实现逻辑非常简单，以下是可直接运行的实现方案，适配你给出的数据集结构。

核心计算逻辑 #

你需要的预测上下边界属于预测区间（区别于对回归均值的置信区间，覆盖单个预测值的波动范围），80%区间对应显著性水平α=0.2，计算步骤如下：

• 拟合模型后，计算训练集残差的无偏标准误
• 对每个待预测点，计算其在特征空间的杠杆值
• 基于t分布（自由度=训练样本数-模型特征数-1）得到对应分位数值
• 边际误差 = t分位数 * 残差标准误 * sqrt(1 + 杠杆值)
• 上下边界 = 预测值 ± 边际误差

代码实现（适配你的数据集） #

首先导入依赖：

import pandas as pd
import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures
from scipy import stats

加载示例数据集，构造时间特征（将日期转为连续序号作为模型输入）：

# 构造示例数据集
df = pd.DataFrame({
    'date': pd.date_range(start='2022-01-01', periods=10, freq='D'),
    'cost': [100,104,107,108,111,117,120,122,128,133]
})
# 构造连续时间特征t，从1开始计数
df['t'] = np.arange(1, len(df)+1)

# 构造未来10天的待预测数据集
future_days = 10
future = pd.DataFrame({
    'date': pd.date_range(start=df['date'].iloc[-1] + pd.Timedelta(days=1), periods=future_days, freq='D'),
    't': np.arange(len(df)+1, len(df)+1+future_days)
})

# 配置80%区间参数
alpha = 0.2

1. 普通线性回归的区间计算 #

# 拟合线性回归模型
X_train = df[['t']].values
y_train = df['cost'].values
X_future = future[['t']].values

lr = LinearRegression()
lr.fit(X_train, y_train)

# 计算点预测值
yhat = lr.predict(X_future)

# 计算80%预测区间
n = len(X_train)
p = X_train.shape[1] # 输入特征数（不含截距）
# 计算残差、残差无偏标准误
y_train_pred = lr.predict(X_train)
resid = y_train - y_train_pred
se = np.sqrt(np.sum(resid**2) / (n - p - 1))
# 计算t分布分位数
t_val = stats.t.ppf(1 - alpha/2, df = n - p -1)
# 补充截距项构造完整特征矩阵，计算叉积逆用于杠杆值计算
X_train_wi = np.hstack([np.ones((n, 1)), X_train])
XtX_inv = np.linalg.inv(X_train_wi.T @ X_train_wi)
# 逐点计算边际误差
margin_err = []
for x in X_future:
    x_wi = np.hstack([[1], x]) # 待预测点补充截距项
    leverage = x_wi @ XtX_inv @ x_wi.T
    me = t_val * se * np.sqrt(1 + leverage)
    margin_err.append(me)
margin_err = np.array(margin_err)

# 汇总结果
lr_result = future.copy()
lr_result['yhat'] = yhat
lr_result['yhat_lower'] = yhat - margin_err
lr_result['yhat_upper'] = yhat + margin_err

2. 多项式回归的区间计算 #

多项式回归仅需先做特征升维，后续区间计算逻辑和普通线性回归完全一致，以2次多项式为例：

# 多项式特征转换，degree参数可根据需求调整
poly = PolynomialFeatures(degree=2, include_bias=False)
X_train_poly = poly.fit_transform(df[['t']].values)
X_future_poly = poly.transform(future[['t']].values)

pr = LinearRegression()
pr.fit(X_train_poly, y_train)
yhat_pr = pr.predict(X_future_poly)

# 计算80%预测区间
n = len(X_train_poly)
p = X_train_poly.shape[1] # 多项式特征数（不含截距）
y_train_pred_pr = pr.predict(X_train_poly)
resid_pr = y_train - y_train_pred_pr
se_pr = np.sqrt(np.sum(resid_pr**2)/(n-p-1))
t_val_pr = stats.t.ppf(1-alpha/2, df = n-p-1)
# 补充截距项构造完整特征矩阵
X_train_poly_wi = np.hstack([np.ones((n,1)), X_train_poly])
XtX_inv_pr = np.linalg.inv(X_train_poly_wi.T @ X_train_poly_wi)

margin_err_pr = []
for x in X_future_poly:
    x_wi = np.hstack([[1], x])
    leverage = x_wi @ XtX_inv_pr @ x_wi.T
    me = t_val_pr * se_pr * np.sqrt(1 + leverage)
    margin_err_pr.append(me)
margin_err_pr = np.array(margin_err_pr)

# 汇总结果
pr_result = future.copy()
pr_result['yhat'] = yhat_pr
pr_result['yhat_lower'] = yhat_pr - margin_err_pr
pr_result['yhat_upper'] = yhat_pr + margin_err_pr

注意事项 #

• 如果需要的是回归均值的80%置信区间，只需要把边际误差公式里的1 + leverage改成leverage即可，计算出的区间会比预测区间窄
• 多项式回归的degree不要设置过高，否则外推预测未来时间点时，区间会快速膨胀，预测结果严重失真
• 该计算方法的前提是模型残差满足独立、同方差、正态分布的假设，如果残差存在明显异方差，需要先对数据做转换或者改用加权最小二乘拟合

内容的提问来源于stack exchange，提问作者user19443045