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

人工智能 偏最小二乘法（PLS）

《人工智能》课程论文

论文题目： 偏最小二乘算法（PLS）回归建模
学生姓名： 张帅帅
学 号： 172341392
专 业： 机械制造及其自动化
所在学院： 机械工程学院

年 月 日
目录
偏最小二乘回归 1
摘要 1
§1偏最小二乘回归原理 1
§2一种更简洁的计算方法 5
§3案例分析 6
致谢 15
附件： 16

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

偏最小二乘回归分析建模的具体步骤如下：
（1）分别提取两变量组的第一对成分，并使之相关性达最大。
（2）假设从两组变量分别提出第一对t₁和u₁,t₁是自变量集,的线性组合：,u₁是因变量集的线性组合：。为了回归分析的需要，要求： 
① t1和u1各自尽可能多地提取所在变量组的变异信息；
② t1和u1的相关程度达到最大。 
由两组变量集的标准化观测数据阵和,可以计算第一对成分的得分向量，记
为和:

第一对成分和的协方差可用第一对成分的得分向量和的内积来计算。故而以上两个要求可化为数学上的条件极值问题：

利用Lagrange乘数法，问题化为求单位向量和,使最大。问题的求解只须通过计算矩阵的特征值和特征向量，且M的最大特征值为,相应的单位特征向量就是所求的解，而可由计算得到。
（3）建立,对的回归及,对的回归。
假定回归模型为：

其中分别是多对一的回归模型中的参数向量，和是残差阵。回归系数向量的最小二乘估计为:

称为模型效应负荷量。
（4）用残差阵和代替和重复以上步骤。
记则残差阵。如果残差阵中元素的绝对值近似为0，则认为用第一个成分建立的回归式精度已满足需要了，可以停止抽取成分。否则用残差阵和代替和重复以上步骤即得：  
分别为第二对成分的权数。而为第二对成分的得分向量。
分别为X,Y的第二对成分的负荷量。这时有

(5)设n×m数据阵的秩为r<=min(n-1,m),则存在r个成分，使得:

把代入，即得p个因变量的偏最小二乘回归方程式：

(6) 交叉有效性检验。
一般情况下，偏最小二乘法并不需要选用存在的r个成分来建立回归
式，而像主成分分析一样，只选用前l个成分（l≤r），即可得到预测能力较好的回归模型。对于建模所需提取的主成分个数l，可以通过交叉有效性检验来确定。  
每次舍去第i个观测，用余下的n-1个观测值按偏最小二乘回归
方法建模，并考虑抽取h个成分后拟合的回归式，然后把舍去的第i个观测点代入所拟合的回归方程式，得到在第i个观测点代入所拟合的预测值。对i=1,2，…n重复以上的验证，即得抽取h个成分时第j个因变量的预测误差平方和为：

的预测误差平方和为：
。
另外，再采用所有的样本点，拟合含h个成分的回归方程。这时，记第i个样本点的预测值为，则可以定义的误差平方和为：

定义Y的误差平方和为：

当达到最小值时，对应的h即为所求的成分个数。通常，总有大于，而则小于。因此，在提取成分时，总希望比值，越小越好；一般可设定限制值为0.05，即当时，增加成分有利于模型精度的提高。或者反过来说，当时，就认为增加新的成分，对减少方程的预测误差无明显的改善作用。
为此，定义交叉有效性为这样，在建模的每一步计算结束前，均进行交叉有效性检验，如果在第h步有则模型达到精度要求，可停止提取成分；若表示第h步提取的成分的边际贡献显著，应继续第h+1步计算。
§2  一种更简洁的计算方法 
上面介绍的算法原则和推导过程的思路在目前的文献中是最为常见的。然而，还有一种更为简洁的计算方法，即直接在矩阵中提取成分 (r≤m)。要求能尽可能多地携带X中的信息，同时，对因变量系统有最大的解释能力。注意，无需在中提取成分得分，这可以使计算过程大为简化，并且对算法结论的解释也更为方便。
偏最小二乘法的简记算法的步骤如下： 
（1） 求矩阵最大特征值所对应的特征向量求得成分计算成分得分向量和残差矩阵其中。
（2） 求矩阵最大特征值所对应的特征向量求得成分计算成分得分向量，和残差矩阵其中。
.
.
.
（r）至第r步，求矩阵最大特征值所对应的特征向量，求得成分计算成分得分向量。  
如果根据交叉有效性，确定共抽取r个成分可以得到一个满意的预测模型,则求在,上的普通最小二乘回归方程为：

把，代入，即得p个因变量的偏最小二乘回归方程式:

这里的满足。
§3  案例分析

本节采用辽宁省给出的关于经济与教育投入历年的数据进行偏最小二乘回归建模。在这个数据系统中被测的样本点，是辽宁省22年的不同教育程度的投资与产业的产出。被测变量分为两组。第一组是自变量包括：L1、L2、L3、L4、K、第二组是因变量包括Y1、Y2、Y3。原始数据见下表1。

表2给出了这8个变量的简单相关系数矩阵。从相关系数矩阵可以看出，自变量与自变量、自变量与因变量、因变量与因变量之间的关系如下：
表2相关系数矩阵
利用如下的MATLAB程序：
%第一步：将样本数据标准化，并求出相关系数矩阵
load jy.txt
mu=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:end);%取因变量数据
e0=data(:,1:n);%取数据标准化的x自变量数值
f0=data(:,n+1:end);%取数据标准化的y因变量数值
num=size(e0,1);%求样本点的个数（也就是说测量的样本多少，本例测量了22年的样本）（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)生成n×n的单位阵
%第二步根据标准化后原始数据矩阵e0和f0计算e0’f0f0’e0的最大特征矩阵所对应的特征向量并计算主元成分ti
%以下计算w，w*和t的主元向量（又称得分向量）
for i=1:n
matrix=e0'*f0*f0'*e0;%根据原始标准化数据矩阵e0，f0，计算w, t p(48)(因变量残差Fi可以不用求p（49）而e0下面循环体依次改变，依次求利用下面循环求wi)
[vec,val]=eig(matrix);%求matrix =e0'*f0*f0'*e0的全部特征值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表示主对角线，k>0表示在主对角线上方，k<0表示在主对角线下方。（在这对角线元素就是特征值λi）
[val,ind]=sort(val,'descend');%降序排列ind表示据单下标换算出全下标
w(:,i)=vec(:,ind(1));%提出最大值对应的特征向量
w_star(:,i)=chg*w(:,i);%计算w*的取值（w*是最大特征值对应的特征向量w*）
t(:,i)=e0*w(:,i);%计算成分t的主元向量（T=E0*W*）p（48）（e0不是固定的在循环体内的）
%第三步建立回归模型，并估计主成分系数pi
pi=e0'*t(:,i)/(t(:,i)'*t(:,i));%计算第i个主成分系数向量pi =pi_i= E0'*ti/(ti'*ti) P(46)-(5-12)
chg=chg*(eye(n)-w(:,i)*pi');%计算w到w*的变换矩阵(w*为用为缩减的自变量数据矩阵X去求新的主元成分ti的对应的权值向量而wi为用为缩减的自变量数据矩阵X的残差矩阵Ei-1去求得ti对应的权值向量eye(n)=I ，I为单位向量)（下次循环用的）p(69) p(51)
%计算数据残差Ei(作为初始矩阵计算下一个成分ti)
e=e0-t(:,i)*pi';%计算残差矩阵
e0=e;%将残差矩阵付给e0,再依次计算下一个主成分（循环计算出所有主成分）
%第四步PLS确定主元r个数采用交叉检验法确定，一般r<m;
%以下计算ss（i）的值即残差的平方和（全部样本的）
%计算回归系数bi(因变量与主成分之间的系数)beta=[t(:,1:i),ones(num,1)]\f0;%求回归方程 的系数（全部样本的因变量与主成分之间的系数不是自变量），数据标准化，没有常数项e=ones(8,1)表示将一个8行1列且元素全为1的矩阵赋值给e
beta(end,:)=[]; %删除回归分析的常数项
cancha=f0- t(:,1:i)*beta;%求矩阵残差矩阵%(e(k)=y(k)-b*t)
ss(i)=sum(sum(cancha.^2));%求误差平方和
%以下求press（i）因变量残差（去掉第j个样本后的预测误差）
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;%求去掉第j个样本后的回归分析的系数，这里带有常数项
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);%p（68）
fprintf(' ');
fprintf('交叉的有效性=%f',Q_h2(i));
r=i
break
end
end
%计算回归系数bi(求Y*关于自变量主元t的回归系数)
beta_z= [t(:,1:r),ones(num,1)]\f0;%求Y*关于自变量主元t的回归系数
beta_z(end,:)=[]; %删除常数项
%第五步根据所求相关回归系数求出自变量Y和X的回归系数，并求出原始回归方程的常数项最后建立回归方程
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);%提出自变量和因变量的标准差
for i=1:m
ch0(i)=mu_y(i)-mu_x./sig_x*sig_y(i)*xishu(:,i); % %计算原始数据的回归方程的常数项
end
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
w1=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%系数矩阵,即未还原原始变量的系数,每一列为一个因变量与自变量的回归方程
计算得提出的成分个数r=4，交叉的有效性=-0.280761；的取值见表3，成分的得分见表4。
表3 的取值
表4 成分的得分
标准化变量关于成分的回归模型如下：

由于成分可以写成原变量的标准化变量的函数，即有

由此可得由成分所建立的偏最小二乘回归模型为：
有关的计算结果见表5.
表5回归系数

所以有：

将标准化变量分别还原成原始变量，则回归方程为：
为了更直观、迅速地观察各个自变量在解释时的边际作用，可以绘制回归系数图，见图1。这个图是针对标准化数据的回归方程的。从回归系数图中可以立刻观察到，这5种自变量那种对产业的影响最大。
为了考察这三个回归方程的模型精度，我们以为坐标值，对所有的样本点绘制预测图。是第k个变量，第i个样本点的预测值。在这个预测图上，如果所有点都能在图的对角线附近均匀分布，则方程的拟合值与原值差异很小，这个方程的拟合效果就是满意的。三种产业的预测图见图2
作图程序如下：
load mydata
ch0=repmat(ch0,num,1);%以ch0的内容堆叠在（numx1）的矩阵ch0
yhat=ch0+x0*xish;%计算Y的预测值
y1max=max(yhat);%求预测值的最大值
y2max=max(y0),%求观测值的最大值
ymax=max([y1max;y2max]);%求预测值和观测值的最大值
cancha=yhat-y0;%计算残差
figure(2)
subplot(2,2,1);%画直线y=x,并画预测图
plot(0:ymax(1),0:ymax(1),yhat(:,1),y0(:,1),'*');
title('第一产业预测')
subplot(2,2,2);
plot(0:ymax(2),0:ymax(2),yhat(:,2),y0(:,2),'O');
title('第二产业预测') ;
subplot(2,1,2);
plot(0:ymax(3),0:ymax(3),yhat(:,3),y0(:,3),'H');
title('第三产业预测')
figure(1)
bar(xishu');
title('回归系数直方图');
% 拟合效果的确定
%所有点都在对角线附近均匀分布，则效果较好

图一

图二

致谢
经过《人工智能》课程的学习，我了解并学习了不同的人工智能算法，包括：深度学习、遗传算法、BP神经网络、蚁群算法等。在林老师的高压教学下，我也自学了PLS偏最小二乘法，并且运用此方法对工程性问题进行了回归建模找到了辽宁省不同产业年产值和不同教育程度的关系，有利于政府加大对其影响因子的投入，从而得到自己想要的结果。在此对本次自学过程中对本次课程作业提供帮助的林老师以及其他同学表示感谢！

附件：
MATLAB程序：
load jy.txt
mu=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:end);%取因变量数据
e0=data(:,1:n);%取数据标准化的x自变量数值
f0=data(:,n+1:end);%取数据标准化的y因变量数值
num=size(e0,1);%求样本点的个数（也就是说测量的样本多少，本例测量了22年的样本）（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)生成n×n的单位阵
%第二步根据标准化后原始数据矩阵e0和f0计算e0’f0f0’e0的最大特征矩阵所对应的特征向量并计算主元成分ti
%以下计算w，w*和t的主元向量（又称得分向量）
for i=1:n
matrix=e0'*f0*f0'*e0;%根据原始标准化数据矩阵e0，f0，计算w, t p(48)(因变量残差Fi可以不用求p（49）而e0下面循环体依次改变，依次求利用下面循环求wi)
[vec,val]=eig(matrix);%求matrix =e0'*f0*f0'*e0的全部特征值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表示主对角线，k>0表示在主对角线上方，k<0表示在主对角线下方。（在这对角线元素就是特征值λi）
[val,ind]=sort(val,'descend');%降序排列ind表示据单下标换算出全下标
w(:,i)=vec(:,ind(1));%提出最大值对应的特征向量
w_star(:,i)=chg*w(:,i);%计算w*的取值（w*是最大特征值对应的特征向量w*）
t(:,i)=e0*w(:,i);%计算成分t的主元向量（T=E0*W*）p（48）（e0不是固定的在循环体内的）
%第三步建立回归模型，并估计主成分系数pi
pi=e0'*t(:,i)/(t(:,i)'*t(:,i));%计算第i个主成分系数向量pi =pi_i= E0'*ti/(ti'*ti) P(46)-(5-12)
chg=chg*(eye(n)-w(:,i)*pi');%计算w到w*的变换矩阵(w*为用为缩减的自变量数据矩阵X去求新的主元成分ti的对应的权值向量而wi为用为缩减的自变量数据矩阵X的残差矩阵Ei-1去求得ti对应的权值向量eye(n)=I ，I为单位向量)（下次循环用的）p(69) p(51)
%计算数据残差Ei(作为初始矩阵计算下一个成分ti)
e=e0-t(:,i)*pi';%计算残差矩阵
e0=e;%将残差矩阵付给e0,再依次计算下一个主成分（循环计算出所有主成分）
%第四步PLS确定主元r个数采用交叉检验法确定，一般r<m;
%以下计算ss（i）的值即残差的平方和（全部样本的）
%计算回归系数bi(因变量与主成分之间的系数)beta=[t(:,1:i),ones(num,1)]\f0;%求回归方程 的系数（全部样本的因变量与主成分之间的系数不是自变量），数据标准化，没有常数项e=ones(8,1)表示将一个8行1列且元素全为1的矩阵赋值给e
beta(end,:)=[]; %删除回归分析的常数项
cancha=f0- t(:,1:i)*beta;%求矩阵残差矩阵%(e(k)=y(k)-b*t)
ss(i)=sum(sum(cancha.^2));%求误差平方和
%以下求press（i）因变量残差（去掉第j个样本后的预测误差）
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;%求去掉第j个样本后的回归分析的系数，这里带有常数项
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);%p（68）
fprintf(' ');
fprintf('交叉的有效性=%f',Q_h2(i));
r=i
break
end
end
%计算回归系数bi(求Y*关于自变量主元t的回归系数)
beta_z= [t(:,1:r),ones(num,1)]\f0;%求Y*关于自变量主元t的回归系数
beta_z(end,:)=[]; %删除常数项
%第五步根据所求相关回归系数求出自变量Y和X的回归系数，并求出原始回归方程的常数项最后建立回归方程
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);%提出自变量和因变量的标准差
for i=1:m
ch0(i)=mu_y(i)-mu_x./sig_x*sig_y(i)*xishu(:,i); % %计算原始数据的回归方程的常数项
end
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
w1=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%系数矩阵,即未还原原始变量的系数,每一列为一个因变量与自变量的回归方程
作图程序如下：
load mydata
ch0=repmat(ch0,num,1);%以ch0的内容堆叠在（numx1）的矩阵ch0
yhat=ch0+x0*xish;%计算Y的预测值
y1max=max(yhat);%求预测值的最大值
y2max=max(y0),%求观测值的最大值
ymax=max([y1max;y2max]);%求预测值和观测值的最大值
cancha=yhat-y0;%计算残差
figure(2)
subplot(2,2,1);%画直线y=x,并画预测图
plot(0:ymax(1),0:ymax(1),yhat(:,1),y0(:,1),'*');
title('第一产业预测')
subplot(2,2,2);
plot(0:ymax(2),0:ymax(2),yhat(:,2),y0(:,2),'O');
title('第二产业预测') ;
subplot(2,1,2);
plot(0:ymax(3),0:ymax(3),yhat(:,3),y0(:,3),'H');
title('第三产业预测')
figure(1)
bar(xishu');
title('回归系数直方图');
% 拟合效果的确定
%所有点都在对角线附近均匀分布，则效果较好

- 20 -