第2讲--numpy--MLR

多元线性回归案例一元线性回归y=a0+a1x实验测定一系列的x和y，求a0、a1min(𝑦𝑖−𝑦𝑖2=min𝑦𝑖−𝑎1𝑥𝑖−𝑎02𝑛𝑖=1𝟏𝒙𝟏𝟏𝒙𝟐……𝟏𝒙𝒏𝒂𝟎𝒂𝟏=𝒚𝟏𝒚𝟐…𝒚𝒏或y=Xa写成矩阵形式一元线性回归计算求解y=Xa已知X和y，求aa=(XTX)-1XTy多元线性回归计算求解𝐲𝟏=𝐚𝟏𝟎+𝐚𝟏𝟏𝐱𝟏+𝐚𝟏𝟐𝐱𝟐+⋯+𝐚𝟏𝐧𝐱𝐧𝐲𝟐=𝐚𝟐𝟎+𝐚𝟐𝟏𝐱𝟏+𝐚𝟐𝟐𝐱𝟐+⋯+𝐚𝟐𝐧𝐱𝐧……𝐲𝐦=𝐚𝐦𝟎+𝐚𝐦𝟏𝐱𝟏+𝐚𝐦𝟐𝐱𝟐+⋯+𝐚𝐦𝐧𝐱𝐧矩阵（一行一个样本，xnk：xn的第k样本）𝟏𝐱𝟏𝟏…𝐱𝐧𝟏𝟏𝐱𝟏𝟐…𝐱𝐧𝟐𝟏………𝟏𝐱𝟏𝐤…𝐱𝐧𝐤𝐚𝟏𝟎𝐚𝟐𝟎…𝐚𝐦𝟎𝐚𝟏𝟏𝐚𝟐𝟏…𝐚𝐦𝟏…………𝐚𝟏𝐧𝐚𝟐𝐧…𝐚𝐦𝐧=𝐲𝟏𝟏𝐲𝟐𝟏…𝐲𝐦𝟏𝐲𝟏𝟐𝐲𝟐𝟐…𝐲𝐦𝟐…………𝐲𝟏𝐤𝐲𝟐𝐤…𝐲𝐦𝐤Y=XA求解：A=(XTX)-1XTYExcel中模拟数据y=a0+a1x自己设定a0和a1的值Excel中，A列名写x，B列名写yA列值自己设定，如用10*rand()B列值按公式给,并加上噪声B列如=0.5+2.2*A2+RAND()/5数据有了，如何求解？a=2.2b=0.5注意：Excel中，每次页面转换时，rand会重新运行，X，Y要同时拷贝矩阵类库numpy的方法有矩阵A，B，将其相乘得到矩阵C：C=A.dot(B)矩阵转置：B=A.T逆：B=np.linalg.inv(A)引用importnumpyasnp矩阵类库numpy的方法B=np.zeros((2,3))Barray([[0.,0.,0.],[0.,0.,0.]])矩阵类库numpy的方法x=np.ones((2,3))xarray([[1.,1.,1.],[1.,1.,1.]])矩阵类库numpy的方法矩阵的+，-运算：C=A+B或D=A-B实矩阵的SVD分解：B=np.linalg.svd(A,full_matrices=False)full_matrices=False一定要写，否则会按复数分解矩阵的伪逆：B=np.linalg.pinv(A)矩阵类库numpy的方法增加新行或列importnumpyasnpa=np.array([[1,2,3.0],[4.0,5,6]])aarray([[1.,2.,3.],[4.,5.,6.]])b=np.ones(3)barray([1.,1.,1.])np.insert(a,0,b,axis=0)array([[1.,1.,1.],[1.,2.,3.],[4.,5.,6.]])矩阵类库numpy的方法两矩阵合并a=np.array([[1,2,3.0],[4.0,5,6]])b=np.ones(2)c=np.c_[b,a]#列合并np.r_[]是行carray([[1.,1.,2.,3.],[1.,4.,5.,6.]])矩阵类库numpy的方法矩阵的行和列数size=A.shapesize[0]是行数，size[1]是列数将一个向量设置为一个对角矩阵s=np.array([1,2,3])只有一个[]号，s是一个向量diag=np.zeros((3,3))diag[:3,:3]=np.diag(s)diag的对角元素是1，2，3从文件加载矩阵Z=np.mafromtxt(c:\\python33\S-093790.txt)Data=Z.dataData就是一个矩阵或者np.savetxt(rd:\333.txt,x)#将x矩阵存储到d:\333.txt文件中y=np.loadtxt(rd:\333.txt)#将存储的矩阵，从文件中读到y矩阵一元线性回归求解设有规律上符合如下方程的一组实验数据y=0.5+2.2x见前面数据模拟多元线性回归设有规律上符合如下方程的一组实验数据y1=1.2x1+0.9x2+3.3x3y2=0.7x1+1.8x2+2.3x3y3=0.4x1+0.7x2+4.9x3自行用随机数模拟数据分析矩阵的公式1.20.93.30.71.82.3系数矩阵Y=XA0.40.74.9A与左侧的矩阵有转置多元线性回归解决的问题系数矩阵Y=XA求解回归系数A，该过程称为建模在A已知时，对于新测Xnew，预报Ynew,称为预报Ynew=XnewA公式推导Y=XAXtY=XtXA(XtX)-1XtY=AE:\pyteach\data\1X.txtE:\pyteach\data\1Y.txt多元线性回归程序importnumpyasnpX=np.mafromtxt(E:\pyteach\data\1X.txt)Y=np.mafromtxt(E:\pyteach\data\1Y.txt)X=X.dataY=Y.dataXt=X.TXtX=Xt.dot(X)XtXinv=np.linalg.inv(XtX)temp=XtXinv.dot(Xt)A=temp.dot(Y)A=A.Tsize=A.shapeforiinrange(size[0]):forjinrange(size[1]):print(A[i][j],\t,end=)print()你专业中，哪些问题属于MLR能否将上述程序编程MLR类MLR类的任务：接受X和Y矩阵根据求系数A的公式，建立模型，既求得A预报：给定新的Xnew，预报Ynew=Xnew*AMLR类importnumpyasnpclassMLR:def__init__(self,X,Y):self.X=Xself.Y=Ydefmodelling(self):Xt=np.transpose(X)XtX=np.dot(Xt,X)XtXinv=np.linalg.inv(XtX)temp=np.dot(XtXinv,Xt)A=np.dot(temp,Y)self.A=AdefoutputCoef(self):A=self.Asize=self.A.shapeforiinrange(size[0]):forjinrange(size[1]):print(A[i][j],\t,end=)print()defpredict(self,X):Y=np.dot(X,self.A)returnY调用MLR类的主函数X=np.mafromtxt(E:\学校教学\python\X.txt)Y=np.mafromtxt(E:\学校教学\python\Y.txt)X=X.dataY=Y.datamlr=MLR(X,Y)mlr.modelling()mlr.outputCoef()Y=mlr.predict(X)print(Y)