偏最小二乘法原理及其matlab代码.pdf-资料库

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

（1）分别提取两变量组的第一对成分，并使之相关性达最大。假设从两组变量分别提出第一对成分为 1t 和 1u ， 1t 是自变量集 X ( 1 = x , , mx T ) 的线性组合： t 1 = xw 11 1 + + xw mm 1 = Xw T 1 ， 1u 是因变量集 Y ( 1 = y , , py T ) 的线性组合： u 1 = yv 11 1 + + yv p 1 p = Yv T 1 。为了回归分析的需要，要求： ① 1t 和 1u 各自尽可能多地提取所在变量组的变异信息； ② 1t 和 1u 的相关程度达到最大。由两组变量集的标准化观测数据阵 0E 和 0F ，可以计算第一对成分的得分向量，记为 1ˆt 和 1ˆu ： 1ˆ t = wE 1 0 = 1ˆ u = vF 10 = n 1 x 11 x ⎡ ⎢ ⎢ ⎢ ⎣ y 11 y n 1 ⎡ ⎢ ⎢ ⎢ ⎣ x m 1 x nm ⎤ ⎥ ⎥ ⎥ ⎦ p y 1 y np ⎤ ⎥ ⎥ ⎥ ⎦ 第一对成分 1t 和 1u 的协方差 1 ut ,(Cov 1 ) ⎤ ⎥ ⎥ ⎥ ⎦ w 11 w m 1 ⎡ ⎢ ⎢ ⎢ ⎣ v 11 v 1 ⎤ ⎥ ⎥ ⎥ ⎦ p ⎡ ⎢ ⎢ ⎢ ⎣ = t 11 n 1 = ⎤ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ t ⎢ ⎣ u ⎡ 11 ⎢ ⎢ u ⎢ ⎣ 可用第一对成分的得分向量 1ˆt 和 1ˆu 的内积 ⎤ ⎥ ⎥ ⎥ ⎦ n 1 来计算。故而以上两个要求可化为数学上的条件极值问题： ˆ,ˆ ut ⎧ < >=< ⎪ 1 1 ⎨ 2 www T = ⎪⎩ 1 1 , vYwE 10 0 vv T = 1 1 1 ,1 >= = xFEw T 1 10 v 1 1 T 0 = 2 ⇒ max 利用Lagrange乘数法，问题化为求单位向量 1w 和 1v ，使 vFEw T Tθ 1 10 1 = 0 ⇒ 最大。问题的求解只须通过计算 mm × 矩阵 EFFEM T= 0 T 0 0 的特征值和特征向量，且 M 的最大特 0 征值为 2 1θ ，相应的单位特征向量就是所求的解 1w ，而 1v 可由 1w 计算得到 v 1 = 1 θ 1 wEF T 0 1 0 。（2）建立 ,1 对 1t 的回归及 y py , ,1 对 1t 的回归。 x mx , 假定回归模型为 -532-

⎧ ⎪ ⎨ ⎪⎩ 0 E F 0 = = ˆ t T α 11 ˆ u T β 11 + + E 1 F 1 其中 m) ααα 1 = 11 ( , , 1 T ， ββ 11 = 1 ( , , T β 1 p ) 分别是多对一的回归模型中的参数向量， 1E 和 1F 是残差阵。回归系数向量 1,βα 的最小二乘估计为 1 ⎧ α ⎪ 1 ⎨ ⎪⎩ β 1 = = ˆ tE T 10 ˆ tF T 10 2 2 ˆ t 1 ˆ t 1 ，称 1,βα 为模型效应负荷量。 1 （3）用残差阵 1E 和 1F 代替 0E 和 0F 重复以上步骤。 ˆ β= F 0 ˆ t T 11 ˆ α= E ˆ t T 11 0 ，记。如果残差阵 1F 中元素的绝对值近似为 0，则认为用第一个成分建立的回归式精度已满足需要了，可以停止抽取成分。否则用残差阵 1E 和 1F 代替 0E 和 0F 重复以上步骤即得：，则残差阵 E 1 F 0 F 1 E ， − = − = 0 0 0 ˆE ˆF w 2 = ( w 21 , , mw 2 T ) ； v 2 = v ( , 21 , pv 2 T ) 分别为第二对成分的权数。而 vF 21 为第二对成分的得分向量。 =β 2 ˆ tF T 1 2 2 ˆ t 2 分别为 YX , 的第二对成分的负荷量。这时有， 2ˆ u = ˆ tET 21 2 ˆ t 2 ， 2ˆ t = wE 2 1 =α 2 ⎧ ⎪ ⎨ ⎪⎩ 0 E F 0 = = ˆ t T T αα 2 11 ˆ t T T ββ 2 11 ˆ t 2 ˆ t 2 + + + + 2 E F 2 （4）设 mn × 数据阵 0E 的秩为 r ≤ min( n − ,1 m ) ，则存在 r 个成分 1 ， t t , 2 rt , , 使得把 t k = 0 = = E F 0 ⎧ ⎪ ⎨ ⎪⎩ 11 xw k + ˆ t T α 11 ˆ t T β 11 + + + + ˆ t T α r r ˆ t T β r r + + r E F r + xw km m （ k ,2,1 = , r ），代入 Y = 11 t β + + rt β r ，即得 p 个因变量的偏最小二乘回归方程式 -533-

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

PRESS h (SS)( −≤−h )05.01( )1 2 = 95.0 2 时，增加成分 ht 有利于模型精度的提高。或者反过来说，当 PRESS h (SS)( >−h 295.0)1 时，就认为增加新的成分 ht ，对减少方程的预测误差无明显的改善作用。为此，定义交叉有效性为 Qh 2 1 −= PRESS h (SS)( h − )1 ，这样，在建模的每一步计算结束前，均进行交叉有效性检验，如果在第 h 步有 2 −

如果根据交叉有效性，确定共抽取 r 个成分 ,1 可以得到一个满意的预测模型， t rt , 则求 0F 在 ,ˆ1 上的普通最小二乘回归方程为 t ˆ, rt F 0 = ˆ t T β 11 + + ˆ t r T β r + F r 把 t k = xw * k 1 1 + + xw * km m （ k ,2,1 = , r ），代入 Y = 1 t β + 1 + rt β r ，即得 p 个因变量的偏最小二乘回归方程式 xa jm 11 xa j + + = y j ，（ j ,2,1 = , m ）。 m h 1 ， ∏− w * h = j 1 = ( wI − α 。 j w h T j ) 这里的 * hw 满足 ˆ t = h wE * h 0 §3 案例分析本节采用兰纳胡德（Linnerud）给出的关于体能训练的数据进行偏最小二乘回归建模。在这个数据系统中被测的样本点，是某健身俱乐部的 20 位中年男子。被测变量分为两组。第一组是身体特征指标 X ，包括：体重、腰围、脉搏。第二组变量是训练结果指标Y ，包括：单杠、弯曲、跳高。原始数据见表 1。 No 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 -536- 表 1 体能训练数据体重(x1) 腰围(x2) 脉搏(x3) 单杠(y1) 弯曲(y2) 跳高(y3) 191 189 193 162 189 182 211 167 176 154 169 166 154 247 193 202 36 37 38 35 35 36 38 34 31 33 34 33 34 46 36 37 5 2 12 12 13 4 8 6 15 17 17 13 14 1 6 12 162 110 101 105 155 101 101 125 200 251 120 210 215 50 70 210 60 60 101 37 58 42 38 40 40 250 38 115 105 50 31 120 50 52 58 62 46 56 56 60 74 56 50 52 64 50 46 62

17 18 19 20 均值标准差 176 157 156 138 178.6 24.6905 37 32 33 33 35.4 3.202 54 52 54 68 56.1 7.2104 4 11 15 2 9.45 5.2863 60 230 225 110 145.55 62.5666 25 80 73 43 70.3 51.2775 表 2 给出了这 6 个变量的简单相关系数矩阵。从相关系数矩阵可以看出，体重与腰围是正相关的；体重、腰围与脉搏负相关；而在单杠、弯曲与跳高之间是正相关的。从两组变量间的关系看，单杠、弯曲和跳高的训练成绩与体重、腰围负相关，与脉搏正相关。表 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); %数据标准化 n=3;m=3; %n 是自变量的个数,m 是因变量的个数 x0=pz(:,1:n);y0=pz(:,n+1:end); e0=data(:,1:n);f0=data(:,n+1:end); 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)); %提出最大特征值对应的特征向量 -537-

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(:,1:i),ones(num,1)]\f0; %求回归方程的系数 beta(end,:)=[]; %删除回归分析的常数项 cancha=f0-t(:,1:i)*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; %求回归分析的系数 beta1(end,:)=[]; %删除回归分析的常数项 cancha=she_f-she_t*beta1; %求残差向量 press_i(j)=sum(cancha.^2); end press(i)=sum(press_i); if i>1 Q_h2(i)=1-press(i)/ss(i-1); else Q_h2(1)=1; end if Q_h2(i)<0.0975 fprintf('提出的成分个数 r=%d',i); r=i; break end end beta_z=[t(:,1:r),ones(num,1)]\f0; %求 Y 关于 t 的回归系数 beta_z(end,:)=[]; %删除常数项 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:end); -538-

资料库

偏最小二乘法原理及其matlab代码.pdf

相关推荐

开发技术

热门标签

最新资料