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

对一元线性回归模型 x = { x 1 , x 2 , ⋯   , x n } x=\{x_1,x_2,\cdots,x_n\} x={x1​,x2​,⋯,xn​}， Y = { Y 1 , Y 2 , ⋯   , Y n } Y=\{Y_1,Y_2,\cdots,Y_n\} Y={Y1​,Y2​,⋯,Yn​}，$Y_i$~$N(a{x}_i+b,{\sigma }^2)，i=1,2,\cdots ,n$，若算得参数$a$，$b$和${\sigma }^2$的估计量$\stackrel{\wedge }{a}$，$\stackrel{\wedge }{b}$和$\stackrel{\wedge }{{\sigma }^2}$。对给定的置信水平$1-\alpha$以及与诸$Y_i,i=1,2,\cdots ,n$独立的随机变量$Y$~$N(ax+b,{\sigma }^2)$的某个取值范围$\Omega$，寻求使得
$$P(Y\in \Omega )\ge 1-\alpha$$
成立的$x$构成的集合其上（下）界的估计量问题，称为控制问题。
例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}$$
要求铸模的抗张强度$Y$的值介于260~340之间，则铝材的硬度应如何控制（置信水平$1-\alpha =0.95$）？就是一个典型的控制问题。
由于$\frac{Y-ax-b}{\sigma }$~$N(0,1)$，用$\stackrel{\wedge }{a}$，$\stackrel{\wedge }{b}$和$\stackrel{\wedge }{{\sigma }^2}$替代$a$，$b$和${\sigma }^2$，$\frac{Y-\stackrel{\wedge }{a}x-\stackrel{\wedge }{b}}{\stackrel{\wedge }{\sigma }}$近似服从$N(0,1)$。对于$\Omega =(y^{\ast },y^{\ast \ast })$的情形，其中$y^{\ast }$和$y^{\ast \ast }$为实数，且满足$y^{\ast \ast }-y^{\ast }>2z_{\alpha /2}\stackrel{\wedge }{\sigma }$，则必有
$$P(y^{\ast }<Y<y^{\ast \ast })=P\left(\frac{y^{\ast }-\stackrel{\wedge }{a}x-\stackrel{\wedge }{b}}{\stackrel{\wedge }{\sigma }}<\frac{Y-\stackrel{\wedge }{a}x-\stackrel{\wedge }{b}}{\stackrel{\wedge }{\sigma }}<\frac{y^{\ast \ast }-\stackrel{\wedge }{a}x-\stackrel{\wedge }{b}}{\stackrel{\wedge }{\sigma }}\right).$$
解不等式$\frac{y^{\ast }-\stackrel{\wedge }{a}x-\stackrel{\wedge }{b}}{\stackrel{\wedge }{\sigma }}\le -z_{\alpha /2}$得$x^{\ast }=\frac{1}{\stackrel{\wedge }{a}}(y^{\ast }-\stackrel{\wedge }{b}+z_{\alpha /2}\stackrel{\wedge }{\sigma })$，解$\frac{y^{\ast \ast }-\stackrel{\wedge }{a}x-\stackrel{\wedge }{b}}{\stackrel{\wedge }{\sigma }}\ge z_{\alpha /2}$得$x^{\ast \ast }=\frac{1}{\stackrel{\wedge }{a}}(y^{\ast \ast }-\stackrel{\wedge }{b}-z_{\alpha /2}\stackrel{\wedge }{\sigma })$。则
$$P(y^{\ast }<Y<y^{\ast \ast })=P\left(\frac{y^{\ast }-\stackrel{\wedge }{a}{x}^{\ast }-\stackrel{\wedge }{b}}{\stackrel{\wedge }{\sigma }}\le -z_{\alpha /2}<\frac{Y-\stackrel{\wedge }{a}x-\stackrel{\wedge }{b}}{\stackrel{\wedge }{\sigma }}<z_{\alpha /2}\le \frac{y^{\ast \ast }-\stackrel{\wedge }{a}{x}^{\ast \ast }-\stackrel{\wedge }{b}}{\stackrel{\wedge }{\sigma }}\right)\ge 1-\alpha .$$
于是，欲使$y^{\ast }<Y<y^{\ast \ast }$，在置信水平$1-\alpha$下，需控制$x\in (x^{\ast },x^{\ast \ast })$（$\stackrel{\wedge }{a}>0$）或$x\in (x^{\ast \ast },x^{\ast })$（$\stackrel{\wedge }{a}<0$）。将上述思想写成如下代码。

from scipy.stats import norm                        #导入norm
def control(a, b, s, y1, y2, alpha):              	#函数定义
    z1,z2=norm.interval(1-alpha)                    #N(0,1)的双侧分位点
    c1=y1-b                                         #y*-b
    c2=y2-b                                         #y**-b
    dy1=z1*s                                        #z1*s
    dy2=z2*s                                        #z2*s
    p1=(c1-dy1)/a                                   #关于y*的端点
    p2=(c2-dy2)/a                                   #关于y**的端点
    if p2<p1:                                       #确定左右端点
        (p1,p2)=(p2,p1)
    return (p1, p2)

程序的第3行计算标准正态分布对应$1-\alpha$的双侧分位点$-z_{\alpha /2}$、$z_{\alpha /2}$，记为z1和z2。第4、5行分别计算$y^{\ast }-\stackrel{\wedge }{b}$和$y^{\ast \ast }-\stackrel{\wedge }{b}$，记为c1和c2。第6、7行分别计算$-z_{\alpha /2}\stackrel{\wedge }{\sigma }$和$z_{\alpha /2}\stackrel{\wedge }{\sigma }$，记为dy1和dy2。第8、9行分别计算$\frac{1}{\stackrel{\wedge }{a}}(y^{\ast }-\stackrel{\wedge }{b}+z_{\alpha /2}\stackrel{\wedge }{\sigma })$和$\frac{1}{\stackrel{\wedge }{a}}(y^{\ast \ast }-\stackrel{\wedge }{b}-z_{\alpha /2}\stackrel{\wedge }{\sigma })$，记为p1和p2。第10~11行的if语句确定控制区间的左、右端点。需要提醒的是，调用函数control前需自行检验$y^{\ast \ast }-y^{\ast }>2z_{\alpha /2}\stackrel{\wedge }{\sigma }$。下列代码完成例1的计算。

import numpy as np                          				#导入numpy
from scipy.stats import linregress          				#导入linregress
alpha=0.05                                  				#设置数据
y1=260
y2=340
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])
n=x.size                                    				#样本容量
x_bar=x.mean()                              				#x数据均值
lxx=((x-x_bar)**2).sum()                    				#lxx
res=linregress(x, y)                        				#调用linregress
a=res.slope                                 				#读取a
b=res.intercept                             				#读取b
s=res.stderr*np.sqrt((n-2)*lxx/n)           				#计算s
print('x in (%.0f, %.0f)'%control(a, b, s, y1, y2, alpha))	#计算控制区间

程序的第3~7行设置原始数据。第9行计算样本容量$n$，第9行计算$x$的数据均值$\overline{x}$记为x_bar。第10行计算$l_{xx}=\sum _{i=1}^n(x_i-\overline{x})$记为lxx。第11行调用函数linregress计算一元回归分析，返回值记为res。第12、13行分别读取$\stackrel{\wedge }{a}$和$\stackrel{\wedge }{b}$，记为a和b。第14行利用res的字段stderr（$=\sqrt{\frac{n\stackrel{\wedge }{{\sigma }^2}}{(n-2)l_{xx}}}$）乘以$\sqrt{\frac{(n-2)l_{xx}}{n}}$，计算$\stackrel{\wedge }{\sigma }$记为s。第15行调用函数contol计算$260<Y<340$的控制区间并输出。运行程序，输出

x in (59, 60)

即若要求铸模的抗张强度$Y$的值介于260~340之间，则铝材的硬度应控制在(59, 60)范围内。
写博不易，敬请支持：
如果阅读本文于您有所获，敬请点赞、评论、收藏，谢谢大家的支持！
代码诚可贵，原理价更高。若为AI学，读正版书好。
返回《导引》