概率统计Python计算：一元线性回归应用——预测

对一元线性回归模型，若算得参数$a$，$b$和${\sigma }^2$的估计量$\stackrel{\wedge }{a}$，$\stackrel{\wedge }{b}$和$\stackrel{\wedge }{{\sigma }^2}$。设$x=x_0$为一指定值，依$E(Y_0)=a{x}_0+b$所得随机变量记为$Y_0$。对置信水平$1-\alpha$，希望寻求统计量$\underline{Y_0}$和$\overline{Y_0}$，使得
$$P(\underline{Y_0}<Y_0<\overline{Y_0})\ge 1-\alpha .$$
这一问题称为预测问题。$(\underline{Y_0},\overline{Y_0})$称为置信水平$1-\alpha$下$Y_0$的预测区间。
由于$\frac{Y_0-\stackrel{\wedge }{a}{x}_0-\stackrel{\wedge }{b}}{\stackrel{\wedge }{\sigma }\sqrt{\frac{n}{n-2}\left[1+\frac{1}{n}+\frac{(x_0-\overline{x}{)}^2}{\sum _{i=1}^n(x_i-\overline{x}{)}^2}\right]}}$~$t(n-2)$，由此可得置信水平$1-\alpha$下$Y_0$的预测区间为
$$\left(\stackrel{\wedge }{a}{x}_0+\stackrel{\wedge }{b}\pm t_{\alpha /2}(n-2)\stackrel{\wedge }{\sigma }\sqrt{\frac{n}{n-2}\left[1+\frac{1}{n}+\frac{(x_0-\overline{x}{)}^2}{\sum _{i=1}^n(x_i-\overline{x}{)}^2}\right]}\right).$$
注意预测区间的增量因子
$$\begin{gathered} \stackrel{\wedge }{\sigma }\sqrt{\frac{n}{n-2}\left[1+\frac{1}{n}+\frac{(x_0-\overline{x}{)}^2}{\sum _{i=1}^n(x_i-\overline{x}{)}^2}\right]}=\stackrel{\wedge }{\sigma }\sqrt{\frac{n+1}{n-2}+\frac{n(x_0-\overline{x}{)}^2}{(n-2)\sum _{i=1}^n(x_i-\overline{x}{)}^2}} \\ =\sqrt{\frac{n+1}{n-2}{\stackrel{\wedge }{\sigma }}^2+\frac{n{\stackrel{\wedge }{\sigma }}^2}{(n-2)\sum _{i=1}^n(x_i-\overline{x}{)}^2}(x_0-\overline{x}{)}^2} \end{gathered}$$
最后的根式内部第2项因子$\frac{n{\stackrel{\wedge }{\sigma }}^2}{(n-2)\sum _{i=1}^n(x_i-\overline{x}{)}^2}$恰为调用linregress函数所得返回值的stderr字段的平方。用linregress函数算得一元回归模型参数$a$，$b$及$\sigma$的无偏估计$\stackrel{\wedge }{a}$，$\stackrel{\wedge }{b}$和$\stackrel{\wedge }{\sigma }$，对给定的$x_0$及置信水平$1-\alpha$，可调用muBounds函数，传递参数mean为$\stackrel{\wedge }{a}{x}_0+\stackrel{\wedge }{b}$，d为$\sqrt{\frac{n+1}{n-2}{\stackrel{\wedge }{\sigma }}^2+\frac{n{\stackrel{\wedge }{\sigma }}^2}{(n-2)\sum _{i=1}^n(x_i-\overline{x}{)}^2}(x_0-\overline{x}{)}^2}$，confidence为$1-\alpha$及df为$n-2$即可求得$Y_0$~$N(a{x}_0+b,{\sigma }^2)$预测区间。
例1设炼铝厂所产铸模的抗张强度与所用铝的硬度有关。设当铝的硬度为$x$时，抗张强度$Y$~$N(ax+b,{\sigma }^2)$，其中$a$，$b$和${\sigma }^2$均未知。对于一系列的$x$值，测得相应的抗张强度如下表
$$\begin{gathered} \text{硬度}x:51,53,60,64,68,70,70,72,83,84 \\ \text{抗张强度}Y:283,293,290,256,288,349,340,354,324,343 \end{gathered}$$
计算对应$x_0=69$，随机变量$Y_0$~$N(a{x}_0+b,{\sigma }^2)$的置信水平为0.95的预测区间。
解： 下列代码完成本例计算。

import numpy as np
from scipy.stats import linregress
x=np.array([51, 53, 60, 64, 68, 70, 70, 72, 83, 84])			#样本数据
y=np.array([283, 293, 290, 286, 288, 349, 340, 354, 324, 343])
alpha=0.05														#显著水平
x0=69															#硬度水平
n=x.size														#样本容量
x_bar=x.mean()													#x的均值
lxx=((x-x_bar)**2).sum()										#x的平方和
res=linregress(x, y)											#计算一元回归
a=res.slope														#a的无偏估计
b=res.intercept													#b的无偏估计
s=res.stderr*np.sqrt((n-2)*lxx/n)								#sigma的最大似然估计
d=np.sqrt((n+1)/(n-2)*s**2+((x0-x_bar)*res.stderr)**2)			#预测区间增量因子
mean=ax0+b														#预测区间中心
confidence=1-alpha												#置信水平
(l, r)=muBounds(mean, d, confidence, df=n-2)					#Y0的预测区间
print('(%.3f, %.3f)'%(l,r))

第3~5行按题面设置各项数据。第6行计算样本容量为n，第7行计算$x$的均值$\overline{x}$为x_bar，第8行计算$\sum _{i=1}^n(x_i-\overline{x}{)}^2$为lxx。第9行调用linregress，计算一元回归分析。第10、11和12行分别读取$a$，$b$及$\sigma$的点估计值$\stackrel{\wedge }{a}$，$\stackrel{\wedge }{b}$和$\stackrel{\wedge }{\sigma }$为a，b和s。第13行计算预测区间的增量因子$\sqrt{\frac{n+1}{n-2}{\stackrel{\wedge }{\sigma }}^2+\frac{n{\stackrel{\wedge }{\sigma }}^2}{(n-2)\sum _{i=1}^n(x_i-\overline{x}{)}^2}(x_0-\overline{x}{)}^2}$为d。第14行计算$\stackrel{\wedge }{a}{x}_0+\stackrel{\wedge }{b}$为mean。第15行计算$1-\alpha$为confidence。第17行调用函数muBounds，计算$Y_0$的预测区间。运行程序，输出

(263.342, 372.259)

即对应硬度$x_0=69$，以0.95的置信水平预测抗张强度$Y_0\in (263.342,372.259)$。
写博不易，敬请支持：
如果阅读本文于您有所获，敬请点赞、评论、收藏，谢谢大家的支持！
代码诚可贵，原理价更高。若为AI学，读正版书好。
返回《导引》