MATLAB——偏最小二乘回归算法

article/2025/10/7 13:02:04

设有q个因变量{y1,y2...yq}和p个自变量{x1,x2...xp}。为了研究因变量和自变量的统计关系，观测n个样本点，构成了自变量与因变量的数据表X= $[x_{1},...,x_{p}]_{n*p}$ 和Y= $[y_{1},...,y_{q}]_{n*q}$ 。部分最小二乘回归分别在X和Y中提取成分 $t_{1}$ 和 $u_{1}$ ,他们分别是x1,...,xp和y1,...,yq的线性组合。提取这两个成分有以下要求：

（1）两个成分尽可能多的携带他们各自数据表的变异信息。

（2）两个成分的相关程度达到最大。

也就是说，他们能够尽可能好地代表各自的数据表，同时自变量成分 $t_{1}$ 对因变量成分 $u_{1}$ 有最强的解释能力。

在第一个成分被成功提取后，分别实施X对 $t_{1}$ 的回归和Y对 $u_{1}$ 的回归。如果回归方程达到满意的精度则终止算法；否则，利用残余信息进行第二轮的成分提取，直到达到一个满意的精度。

·该M文件plsfactor.m是专门提取主函数文件。

%用于提取主函数文件
function [omega,t,pp,XXX,r,YYY]=plsfactor(X0,Y0)
XX=X0'*Y0*Y0'*X0;
[V,D]=eig(XX);
Lamda=max(D);
[MAXLamda,I]=max(Lamda);           %最大特征值对应的特征向量
omega=V(:,I);
%第一主元
t=X0*omega;
pp=X0'*t/(t'*t);
XXX=X0-t*pp';
r=Y0'*t/(t'*t);
YYY=Y0-t*r';

该M文件pls.m是对自变量X和因变量Y进行部分最小二乘回归（偏最小二乘回归）的函数文件

function [beta,VIP]=pls(X,Y)%
[n,p]=size(X);%
[n,q]=size(Y);%
meanX=mean(X);%                           %均值
varX=var(X);         %                          %方差
meanY=mean(Y); %                      
varY=var(Y);          %
%%%数据标准化过程
for i=1:p               %for j=1:n           %X0(j,i)=(X(j,i)-meanX(i))/((varX(i))^0.5);%end
end
for i=1:qfor j=1:nY0(j,i)=(Y(j,i)-meanY(i))/((varY(i))^0.5);%end
end[omega(:,1),t(:,1),pp(:,1),XX(:,:,1),rr(:,1),YY(:,:,1)]=plsfactor(X0,Y0);%
[omega(:,2),t(:,2),pp(:,2),XX(:,:,2),rr(:,2),YY(:,:,2)]=plsfactor(XX(:,:,1),YY(:,:,1));%
PRESShj=0;%
tt0=ones(n-1,2);%for i=1:nYY0(1:(i-1),:)=Y0(1:(i-1),:);%YY0(i:(n-1),:)=Y0((i+1):n,:);%tt0(1:(i-1),:)=t(1:(i-1),:);%tt0(i:(n-1),:)=t((i+1):n,:);%expPRESS(i,:)=(Y0(i,:)-t(i,:)*inv((tt0'*tt0))*tt0'*YY0);%for m=1:qPRESShj=PRESShj+expPRESS(i,m)^2;%end
end
sum1=sum(PRESShj);%
PRESSh=sum(sum1);%for m=1:qfor i=1:nSShj(i,m)=YY(i,m,1)^2;end
end
sum2=sum(SShj);
SSh=sum(sum2);
Q=1-(PRESSh/SSh);
k=3;
%%%循环，提取主元
while Q>0.0975[omega(:,k),t(:,k),pp(:,k),XX(:,:,k),rr(:,k),YY(:,:,k)]=plsfactor(XX(:,:,k-1),YY(:,:,k-1));PRESShj=0;tt00=ones(n-1,k);for i=1:nYY0(1:(i-1),:)=Y0(1:(i-1),:);YY0(i:(n-1),:)=Y0((i+1):n,:);tt00(1:(i-1),:)=t(1:(i-1),:);tt00(i:(n-1),:)=t((i+1):n,:);expPRESS(i,:)=(Y0(i,:)-t(i,:)*((tt00'*tt00)^(-1))*tt00'*YY0);for m=1:qPRESShj=PRESShj+expPRESS(i,m)^2;endendfor m=1:qfor i=1:nSShj(i,m)=YY(i,m,k-1)^2;endendsum2=sum(SShj);SSh=sum(sum2);Q=1-(PRESSh/SSh);if Q>0.0975k=k+1;end
endh=k-1;  %%%提取主元的个数
omegaxing=ones(p,h,q);
%%%还原回归系数
for m=1:qomegaxing(:,1,m)=rr(m,1)*omega(:,1);for i=2:(h)for j=1:(i-1)omegaxingi=(eye(p)-omega(:,j)*pp(:,j)');omegaxingii=eye(p);omegaxingii=omegaxingii*omegaxingi;endomegaxing(:,i,m)=rr(m,i)*omegaxingii*omega(:,i);endbeta(:,m)=sum(omegaxing(:,:,m),2);
end
for i=1:hfor j=1:qrelation(i,j)=sum(prod(corrcoef(t(:,i),Y(:,j))))/2;end
end
%%%%%%%%
Rd=relation.*relation;
RdYt=sum(Rd,2)/q;
Rdtttt=sum(RdYt);
omega22=omega.*omega;
VIP=((p/Rdtttt)*(omega22*RdYt)).^0.5;    %%%计算VIP系数