Python实现均值置信限遇问题，求线性回归代码修正提示

我尝试模拟上图中的均值置信限，但代码存在问题，结果与上图差异较大。我使用置信区间求解斜率和截距，能否有人提供修正提示？谢谢

以下是我的代码：

import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import t
from scipy import stats
x = np.array([3,7,11,15,18,27,29,30,30,31,31, 32,33,33,34,36,36,36,37,38,39,39,39,40,41,42,42,43,44,45,46,47,50])
y = np.array([5,11,21,16,16,28,27,25,35,30,40,32,34,32,34,37,38,34,36,38,37,36,45,39,41,40,44,37,44,46,46,49,51])

mean_x = np.mean(x)
n = len(x)
res = stats.linregress(x, y)
tinv = lambda p, df: abs(t.ppf(p/2, df))
ts = tinv(0.05, n-2)
bins = np.linspace(0,3,54)

plt.plot(x, y, 'o', label='Data Points')
plt.plot(x,res.intercept + res.slope*x, 'c', label='fitted line')
plt.plot(x,(res.intercept+ts*res.intercept_stderr)+(res.slope + ts*res.stderr)*x,'b', label='Upper Limit')
plt.plot(x,(res.intercept-ts*res.intercept_stderr)+(res.slope - ts*res.stderr)*x, 'g' ,label='Lower Limit')
plt.legend()
plt.show()

你当前的核心问题是混淆了「斜率/截距的置信区间」和「均值预测的置信限」：前者是对模型参数的区间估计，会得到两条平行于拟合线的直线；而图中的是针对给定x处预测均值的置信带，是随x变化的曲线（x越接近样本均值，区间越窄），计算逻辑完全不同。

修正思路和代码如下：

1. 先对x排序，避免绘图时线条因x无序出现跳跃
2. 计算每个x对应的拟合均值
3. 针对每个x计算均值置信区间的标准误，再结合t分位数得到上下限

import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import t
from scipy import stats

x = np.array([3,7,11,15,18,27,29,30,30,31,31, 32,33,33,34,36,36,36,37,38,39,39,39,40,41,42,42,43,44,45,46,47,50])
y = np.array([5,11,21,16,16,28,27,25,35,30,40,32,34,32,34,37,38,34,36,38,37,36,45,39,41,40,44,37,44,46,46,49,51])

mean_x = np.mean(x)
n = len(x)
res = stats.linregress(x, y)
tinv = lambda p, df: abs(t.ppf(p/2, df))
ts = tinv(0.05, n-2)

# 对x排序，保证绘图线条连贯
x_sorted = np.sort(x)
# 计算拟合均值
y_hat = res.intercept + res.slope * x_sorted

# 计算均值置信区间的标准误
sum_sq_x = np.sum((x - mean_x)**2)
se_mean = res.stderr * np.sqrt(1/n + (x_sorted - mean_x)**2 / sum_sq_x)
# 得到上下置信限
upper_limit = y_hat + ts * se_mean
lower_limit = y_hat - ts * se_mean

plt.plot(x, y, 'o', label='Data Points')
plt.plot(x_sorted, y_hat, 'c', label='Fitted Line')
plt.plot(x_sorted, upper_limit, 'b', label='Upper Mean Confidence Limit')
plt.plot(x_sorted, lower_limit, 'g', label='Lower Mean Confidence Limit')
plt.legend()
plt.show()

运行这段代码后，你会得到和示例图形态一致的曲线型置信带，这才是均值置信限的正确表现形式。

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