偏最小二乘法回归建模案例

《偏最小二乘法回归建模案例》由会员分享，可在线阅读，更多相关《偏最小二乘法回归建模案例（28页珍藏版）》请在装配图网上搜索。

1、偏最小二乘法回归建模案例人工智能课程论文论文题目： 偏最小二乘算法（PLS）回归建模学生姓名：张帅帅学号： 2专业：机械制造及其自动化所在学院械工程学院年_月_日目录偏最小二乘回归错误!未定义书签。摘要错误!未定义书签。1 偏最小二乘回归原理错误!未定义书签。2 一种更简洁的计算方法错误!未定义书签。3 案例分析 错误!未定义书签。致谢 错误!未定义书签附件： 错误!未定义书签。偏最小二乘回归摘要在实际问题中，经常遇到需要研究两组多重相关变量间的相互依赖关系， 并研究用一组变量（常称为自变量或预测变量）去预测另一组变量（常称为因 变量或响应变量），除了最小二乘准则下的经典多元线性回归分析（ML

2、R）, 提取自变量组主成分的主成分回归分析（PCR ）等方法外，还有近年发展起来的 偏最小二乘（PLS ）回归方法。偏最小二乘回归提供一种多对多线性回归建模的方法，特别当两组变量的 个数很多，且都存在多重相关性，而观测数据的数量（样本量）又较少时，用 偏最小二乘回归建立的模型具有传统的经典回归分析等方法所没有的优点。偏最小二乘回归分析在建模过程中集中了主成分分析，典型相关分析和线 性回归分析方法的特点，因此在分析结果中，除了可以提供一个更为合理的回 归模型外，还可以同时完成一些类似于主成分分析和典型相关分析的研究内 容，提供更丰富、深入的一些信息。本文介绍偏最小二乘回归分析的建模方法；通过例子

3、从预测角度对所建立 的回归模型进行比较。关键词： 主元分析、主元回归、回归建模1 偏最小二乘回归原理考虑 p 个变量 y , y ,.y 与 m 个自变量 x ,x ,.x 的建模问题。偏最小二乘1 2 p 1 2 m回归的基本作法是首先在自变量集中提出第一成分D（D是x，x，x1 2 m的线性组合，且尽可能多地提取原自变量集中的变异信息）；同时在因变量集 中也提取第一成分ui，并要求t占U相关程度达到最大。然后建立因变量y, y ,y与耳的回归，如果回归方程已达到满意的精度，则算法中止。否则1 2 p 继续第二对成分的提取，直到能达到满意的精度为止。若最终对自变量集提取 r个成分t,t.，偏

4、最小二乘回归将通过建立y ,y，y与t,t.的回归1 2 r 1 2 p 1 2 r式，然后再表示为y，y，y与原自变量的回归方程式，即偏最小二乘回归方1 2 p程式。为了方便起见，不妨假定p个因变量y，y，y与m个自变量x，x ,.必 均 1 2 p 1 2 m为标准化变量。因变量组和自变量组的n次标准化观测数据阵分别记为：y11y1 pX 11.X1mF =,E =00yyX .Xn1npn1nm偏最小二乘回归分析建模的具体步骤如下： （1）分别提取两变量组的第一对成分，并使之相关性达最大。（2）假设从两组变量分别提出第一对t和ut是自变量集X =（x，,x，的线1m性组合：t = w x

5、 +. + w x = wtX ,U是因变量集Y二C,., y）的线性组合：1 11 11 m m 1 1 pu二v y +. + v y二VtY。为了回归分析的需要，要求：1 11 11 p p 1 t1和u1各自尽可能多地提取所在变量组的变异信息； t1 和 u1 的相关程度达到最大。由两组变量集的标准化观测数据阵E和F，可以计算第一对成分的得分向量,00记X .11.X1mw11=X .Xwn1t11A 二 E w1 0 1nm1my 11y1 pv11=y yvn1np1pun1U 二 F v1 0 1tn1u11第一对成分t和u的协方差Cov(t,u )可用第一对成分的得分向量i和u

6、 1的内积 1 1 1 1来计算。故而以上两个要求可化为数学上的条件极值问题：= E w , Y v= wtEtF x n maxS 110 10 110 0 1WtW 二 |w |2 二 1, VtV 二 |v |2 二 11 1 1 1 1利用Lagrange乘数法，问题化为求单位向量w和v，使0 = WtEtFV n最大。1 1 1 1 0 0 1问题的求解只须通过计算mxm矩阵M二EtFFtE的特征值和特征向量，且0 0 0 0M的最大特征值为0 2，相应的单位特征向量就是所求的解w，而v可由w计算1 1 1 1v得到二FtE建立i TV，对t1的回归及x x，对的回归。1m假定回归模

7、型为：E t d T + ES 0 1 1 1f u 卩 t + f10 1 1 1其中d =（a，,d）r,卩 卩）分别是多对一的回归模型中的参数向1111m 1111p量，e和f是残差阵。回归系数向量d ,p的最小二乘估计为:1 1 1 1d E Tt /t1S011p:FTt /t101122称d ,p为模型效应负荷量。11(4)用残差阵E和F代替E和F重复以上步骤。1 1 0 0记E = t d t ,F = t p t ,则残差阵E E - E ,F F - F。如果残差 0 1 1 0 1 1 1 0 1 1 0 0F阵1中元素的绝对值近似为0,则认为用第一个成分建立的回归式精度已

8、满足需要了，可以停止抽取成分。否则用残差阵E和F代替E和F重复以上步骤1 1 0 0即得：w =(ww;v =(vv分别为第二对成分的权数。而2 212 m2 212 m11 = Ew ,u2 = Fv为第二对成分的得分向量。1 1 2 2 1 22, PFt t /f221 222=Et t /1 2分别为 X,Y 的第二对成分的负荷量。这时有E t a t +1 a t + E 0 1 1 2 2 2F t 卩 T + t 卩 T + F120 1 1 2 2 2设nxm数据阵E的秩为r=min(n-1,m),则存在r个成分七，使得: 011F t 卩 T +t 卩 T + F011r r

9、r把t w x +. + w x (k 1,2,.,r),代入Y t P +. +1 P ,，即得 p 个因变kk 1 1km m1 1r r量的偏最小二乘回归方程式：y a x +. + a x ,(j 1,2,.,m)jj 1 1jm m(6) 交叉有效性检验。一般情况下，偏最小二乘法并不需要选用存在的r个成分y , y，y来建立回1 2 p归式，而像主成分分析一样，只选用前l个成分(IWr)，即可得到预测能力较好 的回归模型。对于建模所需提取的主成分个数I,可以通过交叉有效性检验来确 定。每次舍去第i个观测(i = 1,2,.n)，用余下的n-1个观测值按偏最小二乘回归方法建模，并考虑抽

10、取h个成分后拟合的回归式，然后把舍去的第i个观测点 代入所拟合的回归方程式，得到在第i个观测点代入所拟合的预测值y (h)。(i)j对i=1,2, .n重复以上的验证，即得抽取h个成分时第j个因变量y (j二i,2,.p)的预测误差平方和为：jPRESSQi)=(y - y (h)2(j = l，2,.p)jij(i ) ji=1Y = (y y )t的预测误差平方和为：1pPRESS(h)= ZPRESS (h)i=1j 。另外，再采用所有的样本点，拟合含h个成分的回归方程。这时，记第i个样 本点的预测值为y (h)，则可以定义y的误差平方和为：ijjSS (h) = f (y - y (h

11、)2ji =1 ij ij定义 Y 的误差平方和为：SS(h)= Zp ss (h)jj=1当PRESS (h)达到最小值时，对应的h即为所求的成分个数。通常，总有 PRESS (h )大于SS (h)，而SS (h)则小于SS (h - 1)。因此，在提取成分 时，总希望比值PRESS (h), SS (h -1)越小越好；一般可设定限制值为， 即当PRESS(h)/SS(h 一 1) = (1 -0.05)2 = 0.952时，增加成分t 有利于h模型精度的提高。或者反过来说，当press (h)/ SS (h- 1)，952时，就认为增加新的成分t ,对减少方程的预测误差无明显的改善作用

12、。h为此，定义交叉有效性为Q2 = 1 - PRESS (h)/SS(h 一1),这样，在建模h的每一步计算结束前，均进行交叉有效性检验，如果在第h步有Q 2 =1 - press (h)/SS (h 一1),则模型达到精度要求，可停止提取成h分；若Q2 1 - 0.952 = 0.0985,表示第h步提取的t成分的边际贡献显 hh著，应继续第 h+1 步计算。2 一种更简洁的计算方法上面介绍的算法原则和推导过程的思路在目前的文献中是最为常见的。然而，还有一种更为简洁的计算方法，即直接在E，E 矩阵中提取成分t，t (r0r -10 rWm)。要求t能尽可能多地携带X中的信息，同时，t对因变量

13、系统F有最h h 0大的解释能力。注意，无需在F中提取成分得分u，这可以使计算过程大为 0h简化，并且对算法结论的解释也更为方便。偏最小二乘法的简记算法的步骤如下：1)2)求矩阵EtF ftE最大特征值所对应的特征向量w,求得成分0 0 0 0 1t = WtX，计算成分得分向量t - E w，和残差矩阵1 1 1 0 1E = E t a T ,其中a = Et t / t1 0 1 1 1 0 12求矩阵EtF ftE最大特征值w,所对应的特征向量求得成分0 0 0 0 1t = WtX，计算成分得分向量，和残差矩阵E - E 12a t,其中1 1 2 1 2 2a Et t / t22

14、 1 2(r)至第r步，求矩阵Et F FtE 最大特征值所对应的特征向量w，求 r 1 0 0 r 1r得成分t - WtX，计算成分得分向量tr - E w。rrrr1 r如果根据交叉有效性，确定共抽取r个成分t，t，t可以得到一个满意的预测1 2 r模型，则求F在11,., t，上的普通最小二乘回归方程为：01rF = t 卩 T + . + t 卩 T + F0 1 1 r r r把t = w* x + . + w* x , (j = 1,2,.r)，代入 Y = t 卩 +. + t 卩， kk 1 1km m1 1r r即得 p 个因变量的偏最小二乘回归方程式:y = a x +.

15、+a x ,(j=12,.m)jj1 1 jm m这里的 w*满足 t = E w*, w* =fi_1 (I - w a T )w。hh0 h hj j hj=13 案例分析世倒 超育脱入对蛭济彩碼溥題他同卫龜繼乌禺析1数据资料及相关槪念為硏究教育投入与产业发展之间存在着具体怎样的相关关系，特收集了如下数据资料一義I辽寺塔1984&005年秋育投入与疑济产出敎摒济斜年協Li-3(-4-K仕Ys15341221561256441951296573961S0.426S.2S9.65S41749552232768939810245074.9328.1115.6E嘟670205 S351741070

16、40161233S392.9357.8154.6:凭了119329394549709680S61124532109.5417.0192.619却31552恥839637753155617J41.94Q2 5246.6翊17633270859SK345932571943951414545.1316.9T.咖1677337&S5S007S59165420L07716K.&540.S353.3150033536571569660M3229033180.S590；34292.:/.J/1245352U44刃務096596254712194.&741.y536.51.993130733615572612

17、630759305120260. S1039.3710.8W41273卿23bO614S6J67X6I98S9P31LAU12594碍g1995142544072&353S7672482439517392.21390.01011.2r.or,162510661137957616449&J9U474,11537,7J 145.9mi231649591666365002554&SS3474.11743,91264.2r濒21264755772439155589256277053L51855.21459.1:f怕2424996465S165644CH264255952D.S2001.51049.4兀

18、门291049834587000722325760719503.42344.4I82L2irnr2971&0271623975胡呢52855043544.42440.62048.1zcn;367472791709233622536991450590.22609.92258.225027阳宛$7KS47355278110S7S5615.S2898*92487.9672611588979222S511771弼了俪7昭43061.62S12.0934214498481590549069162995688143953331733注释表中敎拯摘自辽宁统计年鉴200册.娈量说明衡呈教育投入水平的具体揩标站

19、妬一研究生教育程度（硕士屍博士劳动力數* 人八“一高等教育程度（大学本科及专科）劳动力數平仪 人;厶一中等教育程廈C高中及中专劳动力数*也：2 :厶一初等以下教育程度（小学及文盲）劳动力数厲悻：小；K我育的财政投入 X4i： 7TA.） 经济产出的指标集儿-第一产业（包括林业“牧业、渔业等）产出值 （準隹:忆尢:丘第二产业（包括工业和瑋筑业产出值弊低也壮；為-第三产业（包括济通类的交通运输业、邮电通讯业、商业饮世业、物资供销和 仓储业泾金融.保险业.地质晋查业，房地产、公用事业，居民朋务业.歳游业，咨脚 信息服务业和各类技术服务业，等等）产出值 卑忸：世就 本节采用辽宁省给出的关于经济与教育投

20、入历年的数据进行偏最小二乘回归建模。在这个数据系统中被测的样本点，是辽宁省22 年的不同教育程度的投资与 产业的产出。被测变量分为两组。第一组是自变量包括：L1、L2、L3、L4、K、 第二组是因变量包括Y1、Y2、Y3。原始数据见下表1。義I辽寺吊1984Q005年秋胄投入与轻济产出数擄济斜年協11-sk?M-Kft比作1221561256441951296573961S0.426S.2S9.65S41749552232768959810245074.9328.1115.6670205 S35174107040161233S392.9357.8154.6作了11932939454970968

21、0S61124532109.5417.0192.6茫乳1920315526I5R39637733155617J41.9492.5246,fi17633270859SK345932571943951414545.1316.9倾1677337&S58007559165420L07716K.&540.S353.3如15/M)33530571569660M3229033刘04292一:宀1245352U44575509&S5y96254712J94.&741.y536.51993130733615572612630759305120260. S1039.3710.8萸X127335923bO614S6

22、J67X6I9839P3円12594碍g19?514254407263538767248243951.7392.21390.01011.21962刃服61137957616449&JW474J1537J114S.92316495916663865002554&HS3474.11743,91364.2I濒21264755772439155589256277053L51855.21459.1242t4996465S165644CK264255952D.S2001.51049.42C 口 1291049834587000722325760719503.42344.4182L2ir79S.430fil

23、.62S12.0劄吒934214498481590549069162995688143953331733注释表中敎拯摘自辽宁统计年鉴200册.表2给出了这8个变量的简单相关系数矩阵。从相关系数矩阵可以看出，自变量与自变量、自变量与因变量、因变量与因变量之间的关系如下：表2 相关系数矩阵10.9846630.873672-0.484660.9446660.8642590.889930.889520.98466310.911658-0.494440.9695170.9087820.9243880.9278210.8736720.9116581-0.619620.8944370.8939550.87

24、86190.886966-0.48466-0.49444-0.619621-0.41768-0.44359-0.37609-0.380340.9446660.9695170.894437-0.4176810.9634520.9829440.9871480.8642590.9087820.893955-0.443590.96345210.9831150.9817720.889930.9243880.878619-0.376090.9829440.98311510.9960010.889520.9278210.886966-0.380340.9871480.9817720.9960011利用如下的

25、 MATLAB 程序：%第一步：将样本数据标准化，并求出相关系数矩阵 loadmu=mean(jy); %求平均值sig=std(jy);% 求标准差rr=corrcoef(jy);%求相关系数矩阵data=zscore(jy);% 数据标准化n=5;%n 是自变量个数m=3;%m 是因变量个数x0=jy(:,1:n);% 取自变量数据y0=jy(:, n+1:e nd);%取因变量数据e0=data(:,1:n);%取数据标准化的x自变量数值f0=data(:,n+1:end);%取数据标准化的y因变量数值num=size(e0,1);%求样本点的个数(也就是说测量的样本多少，本例测量了 2

26、2 年的样本)(size(A,n)如果在size函数的输入参数中再添加一项n,并用1或2 为n赋值，则size将返回矩阵的行数或列数。其中r=size(A,1)该语句返回的是矩 阵A的行数，c=size(A,2)该语句返回的是矩阵A的列数)chg=eye(n);%w到w*变换矩阵的初始化eye(n)生成nxn的单位阵 %第二步根据标准化后原始数据矩阵 e0 和 f0 计算 e0f0f0e0 的最大特征矩阵所 对应的特征向量并计算主元成分廿%以下计算w，w*和t的主元向量(又称得分向量)for i=1:nmatrix=eO*fO*fO*eO;%根据原始标准化数据矩阵e0，f0，计算w, t p(

27、48)(因变量 残差Fi可以不用求p (49)而e0下面循环体依次改变，依次求利用下面循环求 wi)vec,val=eig(matrix);%求 matrix =eO*fO*fO*eO 的全部特征值 val 全部特征向量 vec-( V,D=eig(A):求矩阵A的全部特征值，构成对角阵D，并求A的特征向 量构成 V 的全部列向量)val=diag(val);%提出对角线元素X = diag(v,k)当v是一个含有n个元素的向量时， 返回一个n+abs(k)阶方阵X,向量v在矩阵X中的第k个对角线上，k=0表示主 对角线，k0表示在主对角线上方，k0表示在主对角线下方。(在这对角线元 素就是特

28、征值入i)val,ind=sort(val,descend);%降序排列ind表示据单下标换算出全下标w(:,i)=vec(:,i nd(1);%提出最大值对应的特征向量w_star(:,i)=chg*w(:,i);%计算w*的取值(w*是最大特征值对应的特征向量w*) t(:,i)=e0*w(:,i);%计算成分t的主元向量(T=EO*W*) p (48)(e0不是固定的在 循环体内的)%第三步建立回归模型，并估计主成分系数pipi=e0*t(:,i)/(t(:,i)*t(:,i);%计算第 i 个主成分系数向量 pi =pi=EO*ti/(ti*ti) P(46)- (5-12)chg=c

29、hg*(eye(n)-w(:,i)*pi);%计算w到w*的变换矩阵(w*为用为缩减的自变量数 据矩阵X去求新的主元成分廿的对应的权值向量而wi为用为缩减的自变量数据 矩阵X的残差矩阵Ei-1去求得廿对应的权值向量eye(n)=l，I为单位向量)(下 次循环用的) p(69) p(51)%计算数据残差Ei(作为初始矩阵计算下一个成分ti)e=e0-t(:,i)*pi;%计算残差矩阵e0=e;%将残差矩阵付给eO,再依次计算下一个主成分(循环计算出所有主成分)%第四步 PLS 确定主元 r 个数采用交叉检验法确定，一般 r1;Q_h2(i)=1-press(i)/ss(i-1);elseQ_h2

30、(1)=1;endif Q_h2(i)0表示在主对角线上方，k0表示在主对角线下方。(在这对角线元 素就是特征值入i)val,ind=sort(val,descend);%降序排列ind表示据单下标换算出全下标w(:,i)=vec(:,i nd(1);%提出最大值对应的特征向量w_star(:,i)=chg*w(:,i);%计算w*的取值(w*是最大特征值对应的特征向量w*) t(:,i)=eO*w(:,i);%计算成分t的主元向量(T=EO*W*) p (48)(e0不是固定的在 循环体内的)%第三步建立回归模型，并估计主成分系数pipi=e0*t(:,i)/(t(:,i)*t(:,i);%

31、计算第 i 个主成分系数向量 pi =pi=EO*ti/(ti*ti) P(46)- (5-12)chg=chg*(eye(n)-w(:,i)*pi);%计算w到w*的变换矩阵(w*为用为缩减的自变量数 据矩阵X去求新的主元成分廿的对应的权值向量而wi为用为缩减的自变量数据 矩阵X的残差矩阵Ei-1去求得廿对应的权值向量eye(n)=l，I为单位向量)(下 次循环用的) p(69) p(51)%计算数据残差Ei(作为初始矩阵计算下一个成分ti)e=e0-t(:,i)*pi;%计算残差矩阵e0=e;%将残差矩阵付给eO,再依次计算下一个主成分(循环计算出所有主成分)%第四步 PLS 确定主元 r

32、 个数采用交叉检验法确定，一般 r1;Q_h2(i)=1-press(i)/ss(i-1);elseQ_h2(1)=1;endif Q_h2(i)fprintf(提出的成分个数 r=%d,i);%p(68)fprintf( );fpri ntf(交叉的有效性=%f,Q_h2(i);r=ibreakendend%计算回归系数bi(求Y*关于自变量主元t的回归系数)beta_z= t(:,1:r),ones(num,1)fO;%求Y*关于自变量主元t的回归系数 beta_z(end,:)=; %删除常数项%第五步根据所求相关回归系数求出自变量Y和X的回归系数，并求出原始回 归方程的常数项最后建立回

33、归方程xishu= w_star(:,1:r)*beta_z;%求Y*关于X*的回归系数，每一列是一个回归方程 mu_x=mu(1:n);mu_y=mu(n+1:end); %提出自变量和因变量的均值sig_x=sig(1:n );sig_y=sig (n+1:e nd);%提出自变量和因变量的标准差for i=1:mch0(i)=mu_y(i)-mu_x./sig_x*sig_y(i)*xishu(:,i); % %计算原始数据的回归方程的常数项endfor i=1:mxish(:,i)=xishu(:,i)./sig_x*sig_y(i);%计算原始数据回归方程的系数，每一列是一个回 归方

34、程endsol=ch0;xish% %显示回归方程的系数，每一列是一个方程，每一列的第一个数 是常数项,每一列为一个因变量与自变量们的回归方程%此为还原为原始变量后 的方程save mydata x0 y0 num xishu ch0 xishw1=w(:,1)w2=w(:,2)w3=w(:,3)w4=w(:,4)wx1=w_star(:,1)wx2=w_star(:,2)wx3=w_star(:,3)wx4=w_star(:,4)tx1=t(:,1)tx2=t(:,2)tx3=t(:,3)tx4=t(:,4)beta_z %回归系数xishu%系数矩阵，即未还原原始变量的系数，每一列为一个因变量与自变量的回归 方程作图程