Matlab-pls介绍.pdf-资料库

第二十七章偏最小二乘回归分析在实际问题中，经常遇到需要研究两组多重相关变量间的相互依赖关系，并研究用一组变量（常称为自变量或预测变量）去预测另一组变量（常称为因变量或响应变量），除了最小二乘准则下的经典多元线性回归分析（MLR），提取自变量组主成分的主成分回归分析（PCR）等方法外，还有近年发展起来的偏最小二乘（PLS）回归方法。偏最小二乘回归提供一种多对多线性回归建模的方法，特别当两组变量的个数很多，且都存在多重相关性，而观测数据的数量（样本量）又较少时，用偏最小二乘回归建立的模型具有传统的经典回归分析等方法所没有的优点。偏最小二乘回归分析在建模过程中集中了主成分分析，典型相关分析和线性回归分析方法的特点，因此在分析结果中，除了可以提供一个更为合理的回归模型外，还可以同时完成一些类似于主成分分析和典型相关分析的研究内容，提供更丰富、深入的一些信息。本章介绍偏最小二乘回归分析的建模方法；通过例子从预测角度对所建立的回归模型进行比较。 §1 偏最小二乘回归分析 yy 1 L 与 m 个自变量 , , 考虑 p 个因变量 py , 2 2 , , xx mx , 1 L 的建模问题。偏最小二 x 乘回归的基本作法是首先在自变量集中提出第一成分 1t （ 1t 是 ,1 L 的线性组合，且尽可能多地提取原自变量集中的变异信息）；同时在因变量集中也提取第一成分 1u ， y ,1 L 与 1t 的回归，如果回归方并要求 1t 与 1u 相关程度达到最大。然后建立因变量程已达到满意的精度，则算法中止。否则继续第二对成分的提取，直到能达到满意的精度为止。若最终对自变量集提取 r 个成分 rt t t , , 2 , 1 L ，偏最小二乘回归将通过建立 y y t t ,1 L 与原自变量的回归方程式， ,1 L 与 , 2 1 L 的回归式，然后再表示为即偏最小二乘回归方程式。 mx py py py rt , , , , , , 为了方便起见，不妨假定 p 个因变量 y ,1 L 与 m 个自变量化变量。因变量组和自变量组的 n 次标准化观测数据阵分别记为 py , x ,1 L 均为标准 mx , y 11 y 1 x 11 x m 1 ⎡ ⎢ ⎢ ⎢ ⎣ ⎤ ⎥ ⎥ ⎥ ⎦ p = = ， E L L F 0 M y ⎡ ⎢ ⎢ ⎢ ⎣ 偏最小二乘回归分析建模的具体步骤如下：（1）分别提取两变量组的第一对成分，并使之相关性达最大。 X 假设从两组变量分别提出第一对成分为 1t 和 1u ， 1t 是自变量集 y , 线性组合：， 1u 是因变量集 Xw T 1 ⎤ ⎥ ⎥ ⎥ ⎦ M x M y M x L L Y = nm np n 1 n 1 0 ( 1 L= , x py , ) T T ) 的 mx , 的线性组 ( 1 L= t 1 yv 11 1 xw = 11 1 + + + L yv p 1 xw + mm 1 Yv T = 1 u 1 = p L 。为了回归分析的需要，要求：合： ① 1t 和 1u 各自尽可能多地提取所在变量组的变异信息； ② 1t 和 1u 的相关程度达到最大。由两组变量集的标准化观测数据阵 0E 和 0F ，可以计算第一对成分的得分向量，记 -673-

为 1ˆt 和 1ˆu ： 1ˆ t = wE 1 0 = 1ˆ u = vF 10 = n 1 x 11 ⎡ ⎢ ⎢ M x ⎢ ⎣ y ⎡ 11 ⎢ ⎢ ⎢ ⎣ M y n 1 L L L L = t 11 ⎡ ⎢ ⎢ t ⎢ ⎣ u 11 M ⎤ ⎥ ⎥ ⎥ ⎦ n 1 ⎤ ⎥ ⎥ ⎥ ⎦ ⎤ ⎥ ⎥ ⎥ ⎦ x m 1 w 11 M x y 1 ⎡ ⎤ ⎥ ⎢ ⎥ ⎢ M w ⎥ ⎢ ⎦ ⎣ m nm 1 v ⎤ ⎡ ⎤ 11 ⎥ ⎢ ⎥ ⎥ ⎢ ⎥ M M ⎥ ⎢ ⎥ y v ⎦ ⎣ ⎦ np p 1 1 ut ,(Cov 1 p = ⎡ ⎢ ⎢ M u ⎢ ⎣ 可用第一对成分的得分向量 1ˆt 和 1ˆu 的内积 n 1 第一对成分 1t 和 1u 的协方差来计算。故而以上两个要求可化为数学上的条件极值问题： ) ˆ,ˆ ut ⎧ < >=< ⎪ 1 1 ⎨ 2 www T = ⎪⎩ 1 1 , vFwE 10 0 vv T = 1 1 1 ,1 >= = vFEw T 1 10 v 1 1 T 0 = 2 ⇒ max 利用Lagrange乘数法，问题化为求单位向量 1w 和 1v ，使题的求解只须通过计算 mm × 矩阵 EFFEM T= 0 T 0 0 0 征值为 2 1θ ，相应的单位特征向量就是所求的解 1w ，而 1v 可由 1w 计算得到 v 1 = = vFEw T Tθ 1 10 1 最大。问的特征值和特征向量，且 M 的最大特 ⇒ 0 1 θ 1 wEF T 0 1 0 。 x ,1 L 对 1t 的回归。 mx , 0 py （2）建立 y ,1 L 对 1t 的回归及 , 假定回归模型为 ˆ t T α = 11 ˆ u T β = 11 T m) ααα 1 E + 1 F + 1 ββ = 其中 11 1E 和 1F 是残差阵。回归系数向量 2 ⎧ ⎪ ⎨ ⎪⎩ , L E F 0 , ， = 11 ( , 1 1 = ˆ t 1 ˆ t 1 1,βα 为模型效应负荷量。 ˆ tE T 10 ˆ tF T 10 ⎧ α ⎪ 1 ⎨ ⎪⎩ β 1 = 1 ， 2 T , L p ) β 1 ( 1,βα 的最小二乘估计为 1 分别是多对一的回归模型中的参数向量， ˆ α= E 称（3）用残差阵 1E 和 1F 代替 0E 和 0F 重复以上步骤。。如果残差阵 1F F 记， 1 中元素的绝对值近似为 0，则认为用第一个成分建立的回归式精度已满足需要了，可以停止抽取成分。否则用残差阵 1E 和 1F 代替 0E 和 0F 重复以上步骤即得： ˆ β= F 0 ，则残差阵 ˆ t T 11 ˆ t T 11 E 1 F 0 ˆE ˆF E ， − = − = T T 0 0 0 0 ) 2ˆ t = wE 2 1 而， , w w ( = 2 21 2ˆ u = vF 21 ˆ tET 21 2 ˆ t 2 =α 2 , ) v ； mw 2 pv L 2 为第二对成分的得分向量。 L= v 21 ( , , 2 ， =β 2 ˆ tF T 2 1 2 ˆ t 2 -674-

使得 0 L + + ˆ t T α 11 ˆ t T β 11 + E ⎧ = ⎪ ⎨ F ⎪⎩ = 0 xw * + k 1 1 ˆ t T α + r r ˆ t T β + r r k ,2,1 L= （因变量的偏最小二乘回归方程式 xa jm + + L xw * km E F r xa j ，（ L 把 = + + = y j t m m k r j L11 , ) w * hm T L , p ,2,1 L= ） h 1 ， ∏− w * = h ( j 1 = , r ），代入 Y = t β L11 + + rt β r ，即得 p 个分别为 YX , 的第二对成分的负荷量。这时有 0 ⎧ ⎪ ⎨ ⎪⎩ E + F + 2 （4）设 mn × 数据阵 0E 的秩为 r ˆ t T T αα 2 11 ˆ t T T ββ 2 11 E F 0 ˆ t 2 ˆ t 2 = = + + 2 ≤ min( n − ,1 m ) ，则存在 r 个成分 t t , 2 1 L ， rt , , 这里 w * h = ( w * h 1 , 满足 ˆ t = h wE * h 0 wI − α 。 j w h T j ) （5）交叉有效性检验。 , , t t , 2 1 L 来建立回归 l ≤ ），即可得到预测能力较好的回归一般情况下，偏最小二乘法并不需要选用存在的 r 个成分式，而像主成分分析一样，只选用前l 个成分（ r 模型。对于建模所需提取的成分个数l ，可以通过交叉有效性检验来确定。每次舍去第i 个观测（），对余下的 1−n 个观测值用偏最小二乘回归方法建模，并考虑抽取 h 个成分后拟合的回归式，然后把舍去的第i 个观测点代入所拟。对 ,1=i 合的回归方程式，得到 ,2 L 重复以上的验证，即得抽取 h 个成分时第 j 个因变量 p ) 的预测误差平方和为在第i 个观测点上的预测值 ˆ )( y j h )( i ,2,1 , L= ,2,1 L= ,2,1 L= n, y j y j rt n p ) ( ( i j j , , PRESS h )( = j n ∑ i 1 = ( y ij − ˆ y i )( j ( h )) 2 （ j ,2,1 L= , p ） ( 1 L= y , , py T ) 的预测误差平方和为 Y PRESS h )( = PRESS j h )( 。 p ∑ i 1 = 另外，再采用所有的样本点，拟合含 h 个成分的回归方程。这时，记第i 个样本点的预测值为 )(ˆ hyij ，则可以定义 jy 的误差平方和为 h )(SS j = ( y ij − (ˆ hy ij 2)) n ∑ i 1 = 定义Y 的误差平方和为 ∑ h )(SS = p j h )(SS j 1 = 当 PRESS h 大于值 PRESS h 达到最小值时，对应的 h 即为所求的成分个数。通常，总有 (SS −h 。因此，在提取成分时，总希望比 )( PRESS )( )(SS h ，而 −h 越小越好；一般可设定限制值为 0.05，即当 )(SS h 则小于 h (SS)( )1 )1 -675-

PRESS h (SS)( −≤−h )05.01( )1 2 = 95.0 2 时，增加成分 ht 有利于模型精度的提高。或者反过来说，当时，就认为增加新的成分 ht ，对减少方程的预测误差无明显的改善作用。 295.0)1 h (SS)( PRESS >−h 2 = .0 .0 − 2 Qh 2 1 −= 0975 0985 2 ≥hQ PRESS 为此，定义交叉有效性为 h −

果指标Y ，包括：单杠、弯曲、跳高。原始数据见表 1。表 2 给出了这 6 个变量的简单相关系数矩阵。从相关系数矩阵可以看出，体重与腰围是正相关的；体重、腰围与脉搏负相关；而在单杠、弯曲与跳高之间是正相关的。从两组变量间的关系看，单杠、弯曲和跳高的训练成绩与体重、腰围负相关，与脉搏正相关。 No 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 均值标准差体重(x1) 腰围(x2) 脉搏(x3) 单杠(y1) 弯曲(y2) 跳高(y3) 表 1 体能训练数据 191 189 193 162 189 182 211 167 176 154 169 166 154 247 193 202 176 157 156 138 178.6 24.6905 36 37 38 35 35 36 38 34 31 33 34 33 34 46 36 37 37 32 33 33 35.4 3.202 50 52 58 62 46 56 56 60 74 56 50 52 64 50 46 62 54 52 54 68 56.1 7.2104 5 2 12 12 13 4 8 6 15 17 17 13 14 1 6 12 4 11 15 2 9.45 5.2863 162 110 101 105 155 101 101 125 200 251 120 210 215 50 70 210 60 230 225 110 145.55 62.5666 60 60 101 37 58 42 38 40 40 250 38 115 105 50 31 120 25 80 73 43 70.3 51.2775 表 2 相关系数矩阵 1 0.8702 -0.3658 -0.3897 -0.4931 -0.2263 0.8702 1 -0.3529 -0.5522 -0.6456 -0.1915 -0.3658 -0.3529 1 0.1506 0.225 0.0349 -0.3897 -0.5522 0.1506 1 0.6957 0.4958 -0.4931 -0.6456 0.225 0.6957 1 0.6692 -0.2263 -0.1915 0.0349 0.4958 0.6692 1 利用如下的 MATLAB 程序： clc,clear load pz.txt %原始数据存放在纯文本文件 pz.txt 中 mu=mean(pz);sig=std(pz); %求均值和标准差 rr=corrcoef(pz); %求相关系数矩阵 data=zscore(pz); %数据标准化,变量记做 X*和 Y* n=3;m=3; %n 是自变量的个数,m 是因变量的个数 x0=pz(:,1:n);y0=pz(:,n+1:end); %原始的自变量和因变量数据 e0=data(:,1:n);f0=data(:,n+1:end); %标准化后的自变量和因变量数据 -677-

num=size(e0,1);%求样本点的个数 chg=eye(n); %w 到 w*变换矩阵的初始化 for i=1:n %以下计算 w，w*和 t 的得分向量， matrix=e0'*f0*f0'*e0; [vec,val]=eig(matrix); %求特征值和特征向量 val=diag(val); %提出对角线元素，即提出特征值 [val,ind]=sort(val,'descend'); w(:,i)=vec(:,ind(1)); %提出最大特征值对应的特征向量 w_star(:,i)=chg*w(:,i); %计算 w*的取值 t(:,i)=e0*w(:,i); %计算成分 ti 的得分 alpha=e0'*t(:,i)/(t(:,i)'*t(:,i)); %计算 alpha_i chg=chg*(eye(n)-w(:,i)*alpha'); %计算 w 到 w*的变换矩阵 e=e0-t(:,i)*alpha'; %计算残差矩阵 e0=e; %以下计算 ss(i)的值 beta=t\f0; %求回归方程的系数，数据标准化，没有常数项 cancha=f0-t*beta; %求残差矩阵 ss(i)=sum(sum(cancha.^2)); %求误差平方和 %以下计算 press(i) for j=1:num t1=t(:,1:i);f1=f0; she_t=t1(j,:);she_f=f1(j,:); %把舍去的第 j 个样本点保存起来 t1(j,:)=[];f1(j,:)=[]; %删除第 j 个观测值 beta1=[t1,ones(num-1,1)]\f1; %求回归分析的系数,这里带有常数项 cancha=she_f-she_t*beta1(1:end-1,:)-beta1(end,:); %求残差向量 press_i(j)=sum(cancha.^2); %求误差平方和 end press(i)=sum(press_i); Q_h2(1)=1; if i>1, Q_h2(i)=1-press(i)/ss(i-1); end if Q_h2(i)<0.0975 fprintf('提出的成分个数 r=%d',i); break end end beta_z=t\f0; %求 Y*关于 t 的回归系数 xishu=w_star*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:end); %提出自变量和因变量的标准差 ch0=mu_y-(mu_x./sig_x*xishu).*sig_y; %计算原始数据回归方程的常数项 for i=1:m xish(:,i)=xishu(:,i)./sig_x'*sig_y(i); %计算原始数据回归方程的系数 end sol=[ch0;xish] %显示回归方程的系数，每一列是一个方程，每一列的第一个数是常数项 save mydata x0 y0 num xishu ch0 xish -678-

1,t t 即可，交叉有效性计算得只要提出两个成分 2 见表 3，成分 ht 的得分 htˆ 见表 4。 −=Q 2 2 0.3263 。 hw 与 * hw 的取值自变量 x1 x2 x3 表 3 hw 的取值 * hw 与 w2 -0.4688 0.5680 0.6765 w1 -0.5899 -0.7713 0.2389 * 1w -0.5899 -0.7713 0.2389 * 2w -0.3679 0.6999 0.6356 表 4 成分 ht 的得分 htˆ 6 4 5 1 2 3 7 11 12 13 15 0.76 -0.749 0.5211 0.2034 -0.3929 -1.4037 -0.8232 -4.3898 -0.4867 -0.9738 -1.1328 -0.2291 1.1867 0.7433 0.3645 0.68 14 0.0767 17 0.0717 16 -0.7007 -0.6983 0.757 -0.5914 -0.1667 0.5212 -0.6429 -0.7697 -0.9074 0.6884 No ˆt 1 ˆt 2 No ˆt 1 ˆt 2 标准化变量 ky~ 关于成分 2 1,t t 的回归模型如下 3,2,1=k 由于成分 ht 可以写成原变量的标准化变量 jx~ 的函数，即有 ~ ~ xwxw * + h h 1 1 2 t 所建立的偏最小二乘回归模型为 1,t 由此可得由成分 2 ~ ~ )~ ~ ~ ~ y xw xwxwr xwxwr ( ( * * * * = + + + + + k k k 1 3 31 2 1 2 1 11 2 22 21 ~) ~) xwr wr xwr wr wr ( ( ( * * * * = + + + k k k k k 11 1 21 1 1 2 2 31 1 12 r r r ( , , ) 的计算结果见表 5。 h h h 1 表 5 回归系数 hr ~ xw * h 3 3 2,1=h * 12 * 22 r = h tr k 11 tr k 2 有关 ~ y k ， 2 + 2 t h = 2 3 ， * 2 + = + k r1 r2 −= −= −= 所以，有 ~ y .0 1 ~ y .0 2 ~ y .0 3 将标准化变量程为 y 1 y = = .47 .612 2 ~ x 0778 .0 − 1 ~ x 1385 .0 − 1 ~ x 0604 .0 − 1 (~,~ =kxy k k ~ x 4989 2 ~ x 5244 2 ~ x 1559 2 )3,2,1 分别还原成原始变量 ~ x 1322 3 ~ x 0854 3 ~ x 0073 3 .0 .0 .0 − − − x 0197 .0 0167 − 1 x 5671 .0 3509 − 1 .0 8237 − .10 − x 2 2477 − x 2 .0 − 0969 .0 x 3 7412 x 3 1 0.3416 -0.3364 2 0.4161 -0.2908 3 0.1430 -0.0652 8 9 10 0.7436 1.7151 1.1626 0.2106 0.6549 -0.1668 18 19 20 1.1993 1.0485 1.9424 -0.7827 -0.3729 1.1294 )~ xw * 32 3 ~) xwr + k 3 2 * 32 =kxy k ( , k )3,2,1 ，则回归方 -679-

= − x 1 − y 3 .2 4969 x .183 .0 1253 9849 x 3 )3,2,1 时的边际作用，可以绘 0518 .0 − 2 ( =kyk 为了更直观、迅速地观察各个自变量在解释制回归系数图，见图 1。这个图是针对标准化数据的回归方程的。从回归系数图中可以立刻观察到，腰围变量在解释三个回归方程时起到了极为重要的作用。然而，与单杠及弯曲相比，跳高成绩的回归方程显然不够理想，三个自变量对它的解释能力均很低。 ,ˆ( ik y y 为了考察这三个回归方程的模型精度，我们以为坐标值，对所有的样本点绘制预测图。 ikyˆ 是第 k 个变量，第i 个样本点 ) ( iky 的预测值。在这个预测图上，如果所有点都能在图的对角线附近均匀分布，则方程的拟合值与原值差异很小，这个方程的拟合效果就是满意的。体能训练的预测图见图 2。 ) ik 0 -0.1 -0.2 -0.3 -0.4 -0.5 -0.6 -0.7 单杠 1 弯曲 2 跳高 3 图 1 回归系数的直方图 300 200 100 0 0 100 200 300 弯曲成绩预测图 20 15 10 5 0 0 250 200 150 100 50 0 0 5 10 15 20 单杠成绩预测图 100 200 300 跳高成绩预测图图 2 体能训练预测图画直方图的 MATLAB 程序为：bar(xishu') 画体能训练的预测图的 MATLAB 程序如下： load mydata ch0=repmat(ch0,num,1); yhat=ch0+x0*xish; %计算 Y 的预测值 y1max=max(yhat); %求预测值的最大值 y2max=max(y0); %求观测值的最大值 ymax=max([y1max;y2max]) %求预测值和观测值的最大值 cancha=yhat-y0; %计算残差 subplot(2,2,1) -680-

资料库

Matlab-pls介绍.pdf

相关推荐

课程资源

热门标签

最新资料