MATLAB——偏最小二乘迴歸演算法
阿新 • • 發佈:2019-01-02
設有q個因變數{y1,y2...yq}和p個自變數{x1,x2...xp}。為了研究因變數和自變數的統計關係,觀測n個樣本點,構成了自變數與因變數的資料表X= 和Y= 。部分最小二乘迴歸分別在X和Y中提取成分和,他們分別是x1,...,xp和y1,...,yq的線性組合。提取這兩個成分有以下要求:
(1)兩個成分儘可能多的攜帶他們各自資料表的變異資訊。
(2)兩個成分的相關程度達到最大。
也就是說,他們能夠儘可能好地代表各自的資料表,同時自變數成分對因變數成分有最強的解釋能力。
在第一個成分被成功提取後,分別實施X對的迴歸和Y對的迴歸。如果迴歸方程達到滿意的精度則終止演算法;否則,利用殘餘資訊進行第二輪的成分提取,直到達到一個滿意的精度。
·該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:q for j=1:n Y0(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:n YY0(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:q PRESShj=PRESShj+expPRESS(i,m)^2;% end end sum1=sum(PRESShj);% PRESSh=sum(sum1);% for m=1:q for i=1:n SShj(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:n YY0(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:q PRESShj=PRESShj+expPRESS(i,m)^2; end end for m=1:q for i=1:n SShj(i,m)=YY(i,m,k-1)^2; end end sum2=sum(SShj); SSh=sum(sum2); Q=1-(PRESSh/SSh); if Q>0.0975 k=k+1; end end h=k-1; %%%提取主元的個數 omegaxing=ones(p,h,q); %%%還原迴歸係數 for m=1:q omegaxing(:,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; end omegaxing(:,i,m)=rr(m,i)*omegaxingii*omega(:,i); end beta(:,m)=sum(omegaxing(:,:,m),2); end for i=1:h for j=1:q relation(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係數