chengxu-plsr2.txt

偏最小二乘程序 标准化 求矩阵的最大特征值及对应特征向量

function [c,m,v]=norm1(C) 
%对数据进行标准化处理 
[n,s]=size(C); 
for i=1:n 
for j=1:s 
c(i,j)=(C(i,j)-mean(C(:,j)))/sqrt(cov(C(:,j))); 
end 
end 
m=mean(C); 
for j=1:s 
v(1,j)=sqrt(cov(C(:,j))); 
end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function w=maxdet(A) 
%求矩阵的最大特征值 
[v,d]=eig(A); 
[n,p]=size(d); 
d1=d*ones(p,1); 
d2=max(d1); 
i=find(d1==d2); 
w=v(:,i); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [t,q,w,wh,f0,FF]=fun717(px,py,C) 
% px自变量的输入个数 
% py因变量的输入个数。 
% C输入的自变量和因变量组成的矩阵 
% t提取的主成分 
% q为回归系数。 
% w最大特征值所对应的特征向量。 
% wh处理后的特征向量 
% f0回归的标准化的方程系数 
% FF原始变量的回归方程的系数 
c=norm1(C); %norm1为标准化函数 
%y=c(:,px+1:px+py); %截取标准化的因变量 
E0=c(:,1:px); %截取标准化的自变量
F0=c(:,px+1:px+py); %截取标准化的因变量
A=E0'*F0*F0'*E0; %关联矩阵
w(:,1)=maxdet(A); % maxdet为求取最大特征值的函数
t(:,1)=E0*w(:,1); %第一主成分 
p(:,1:px)=(E0'*t(:,1)/(t(:,1)'*t(:,1)))'; % p参
E(:,1:px)=E0-t(:,1)*(E0'*t(:,1)/(t(:,1)'*t(:,1)))'; % 残差矩阵
%进入循环
for i=0:px-2 
B(:,px*i+1:px*i+px)=E(:,px*i+1:px*i+px)'*F0*F0'*E(:,px*i+1:px*i+px); 
w(:,i+2)=maxdet(B(:,px*i+1:px*i+px));  
t(:,i+2)=E(:,px*i+1:px*i+px)*w(:,i+2); 
p(:,px*i+px+1:px*i+2*px)=(E(:,px*i+1:px*i+px)'*t(:,i+2)/(t(:,i+2)'*t(:,i+2)))'; 
E(:,px*i+px+1:px*i+2*px)=E(:,px*i+1:px*i+px)-t(:,i+2)*(E(:,px*i+1:px*i+px)'*t(:,i+2)/(t(:,i+2)'*t(:,i+2)))'; 
end 

for s=1:px 
q(:,s)=p(1,px*(s-1)+1:px*s)'; 
end 

[n,d]=size(q); 

for h=1:px 
iw=eye(d); 
for j=1:h-1 
iw=iw*(eye(d)-w(:,j)*q(:,j)'); 
end 
wh(:,h)=iw*w(:,h); 
end

for j=1:py 
zr(j,:)=(regress(F0(:,j),t))';
end 

for j=1:px 
for i=1:py 
w1=wh(:,1:j); 
zr1=(zr(i,1:j))'; 
f0(i,:,j)=(w1*zr1)'; %生成标准化变量的方程的系数矩阵 
end 
[c,m,v]=norm1(C); 
ccxx=ones(py,1)*m(1,1:px); 
ccy=(v(1,px+1:px+py))'*ones(1,px); 
ccx=ones(py,1)*(v(1,1:px)); 
ff=ccy.*f0(:,:,j)./ccx; 
fff=-(sum((ccy.*ccxx.*f0(:,:,j)./ccx)')-m(1,px+1:px+py))'; 
FF(:,:,j)=[fff,ff]; %生成原始变量方程的常数项和系数矩阵 
end 


实验结果：
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[t,q,w,wh,f0,FF]=fun717(3,3,C)

t =

   -0.6429   -0.5914   -0.1305
   -0.7697   -0.1667    0.1339
   -0.9074    0.5212    0.0477
    0.6884    0.6800    0.3458
   -0.4867   -1.1328   -0.1818
   -0.2291    0.0717    0.0247
   -1.4037    0.0767   -0.5718
    0.7436    0.2106   -0.0322
    1.7151    0.6549   -1.5572
    1.1626   -0.1668    0.3329
    0.3645   -0.7007    0.2011
    0.7433   -0.6983    0.0023
    1.1867    0.7570    0.3362
   -4.3898    0.7600    0.2554
   -0.8232   -0.9738   -0.0826
   -0.7490    0.5211   -0.6672
   -0.3929    0.2034    0.5642
    1.1993   -0.7827    0.0924
    1.0485   -0.3729    0.3190
    1.9424    1.1294    0.5677


q =

   -0.6659    0.0198   -0.6575
   -0.6760    0.3547    0.2871
    0.3589    1.1942   -0.6967


w =

   -0.5899   -0.4688   -0.6575
   -0.7713    0.5680    0.2871
    0.2389    0.6765   -0.6967


wh =

   -0.5899   -0.3679   -0.9346
   -0.7713    0.6999    0.8023
    0.2389    0.6356   -0.2228


f0(:,:,1) =

   -0.2015   -0.2635    0.0816
   -0.2454   -0.3209    0.0994
   -0.0843   -0.1103    0.0342


f0(:,:,2) =

   -0.0778   -0.4989   -0.1322
   -0.1385   -0.5244   -0.0854                       %取两个成分的标准数据回归系数
   -0.0604   -0.1559   -0.0073


f0(:,:,3) =

    0.3683   -0.8818   -0.0258
    0.2872   -0.8898    0.0161
   -0.2590    0.0146   -0.0546


FF(:,:,1) =

   29.2002   -0.0431   -0.4350    0.0598
  430.2511   -0.6220   -6.2712    0.8625
  150.4807   -0.1752   -1.7662    0.2429


FF(:,:,2) =

   47.0197   -0.0167   -0.8237   -0.0969
  612.5671   -0.3509  -10.2477   -0.7412            %取两个成分的原始数据回归系数
  183.9849   -0.1253   -2.4969   -0.0518


FF(:,:,3) =

   47.9684    0.0788   -1.4558   -0.0190
  623.2817    0.7277  -17.3872    0.1393
  179.8868   -0.5379    0.2338   -0.3886


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%主成分分析程序
px=3;
py=1;
c=norm1(A);
E0=c(:,1:px);
F0=c(:,px+py);
f0=pinv(E0'*E0)*E0'*F0;
disp(['y1=' num2str(f0(1)) 'x1+(' num2str(f0(2)) ')x2+(' num2str(f0(3)) ')x3'])
[c,m,v]=norm1(A); 
%待改进
for j=1:py 
ccxx=ones(py,1)*m(1,1:px); 
ccy=(v(1,px+1:px+py))'*ones(1,px); 
ccx=ones(py,1)*(v(1,1:px)); 
ff=ccy.*f0(:,:,j)./ccx; 
fff=-(sum((ccy.*ccxx.*f0(:,:,j)./ccx)')-m(1,px+1:px+py))'; 
FF(:,:,j)=[fff,ff];
end 
FF
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function [r,Rdyt,RdYt,RdYtt,Rdytt,VIP]=fun8y(px,py,c) 
X=c(:,1:px); 
Y=c(:,px+1:px+py); 
x=norm1(X); 
y=norm1(Y); 
[t,q,w]=fun717(px,py,[X,Y]); 
r1=corrcoef([y,t]); 
r=r1(py+1:px+py,1:py)'; 
Rdyt=r.^2; 
RdYt=mean(Rdyt) 
for m=1:px 
RdYtt(1,m)=sum(RdYt(1,1:m)'); 
end 
for j=1:py 
for m=1:py 
Rdytt(j,m)=sum(Rdyt(j,1:m)'); 
end 
end 
for j=1:px 
for m=1:px 
Rd(j,m)=RdYt(1,1:m)*((w(j,1:m).^2)'); 
end 
end 
for j=1:px 
VIP(j,:)=sqrt((px*ones(1,px)./RdYtt).*Rd(j,:)); 
end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function [r,Rdxt,RdXt,RdXtt,Rdxtt]=fun8x(px,py,c) 
X=c(:,1:px); 
Y=c(:,px+1:px+py); 
x=norm1(X); 
y=norm1(Y); 
[t,q,w]=fun717(px,py,[X,Y]); 
r1=corrcoef([x,t]); 
r=r1(px+1:px+px,1:px)'; 
Rdxt=r.^2; 
RdXt=mean(Rdxt); 
for m=1:px 
RdXtt(1,m)=sum(RdXt(1,1:m)'); 
end 
for j=1:px 
for m=1:px 
Rdxtt(j,m)=sum(Rdxt(j,1:m)'); 
end 
end 
% for j=1:px 
% for m=1:px 
% Rd(j,m)=RdXt(1,1:m)*((w(j,1:m).^2)'); 
% end 
% end 
% for j=1:px 
% VIP(j,:)=sqrt((px*ones(1,px)./RdYtt).*Rd(j,:)); 
% end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function [t,u]=TU(px,py,C) 
%t提取的自变量的主成分 
%u 提取的因变量的主成分 
c=norm1(C); 
y=c(:,px+1:px+py); 
E0=c(:,1:px); 
F0=c(:,px+1:px+py); 
A=E0'*F0*F0'*E0; 
w(:,1)=maxdet(A); 
t(:,1)=E0*w(:,1); 
B=F0'*E0*E0'*F0; 
cc(:,1)=maxdet(B); 
u(:,1)=F0*cc(:,1); 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function drew(px,py,c) 
X=c(:,1:px); 
Y=c(:,px+1:px+py); 
[line,l]=size(Y); 
[t,q,w,wh,f0,FF]=fun717(px,py,c); 
YY=X*FF(:,2:px+1,3)'+ones(line,1)*FF(:,1,3)'; 

subplot(1,1,1,1) 
bar(f0(:,:,3)) 
title(' 直方图') 
legend('SG','TZBFB','FHL','JK','HPZD','JPZD','TZ','ZG','GPK') 
grid on 
plot(YY(:,4),Y(:,4),'+'); 
lsline 

for i=1:py 
v=mod(i,4); 
d=(i-v)/4; 
subplot(2,2,v,d+1) 
plot(YY(:,i),Y(:,i),'*'); 
lsline 
end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function [ Qhj,Qh,prey]=crossval7(px,py,c) 
%px是自变量的个数； 
%py是因量 
PRESShj=zeros(px,py); 
X=c(:,1:px); 
Y=c(:,px+1:px+py); 
x=norm1(X); 
y=norm1(Y); 
[line,row]=size(x); 
for h=1:px 
for j=1:line 
newx=X; 
newy=Y; 
newx(j,:)=[]; 
newy(j,:)=[]; 
[t,p0,w,wh,f0,FF]=fun717(px,py,[newx,newy]); 
prey(j,:,h)=X(j,:)*FF(:,2:px+1,h)'+FF(:,1,h)'; 
end 
PRESShj(h,:)=sum((Y-prey(:,:,h)).^2); 
end 
PRESSh=(sum(PRESShj'))'; 

for h=1:px 
[t1,p0,w,wh,f0,FF]=fun717(px,py,c); 
prey2(:,:,h)=X(:,:)*FF(:,2:px+1,h)'+ones(line,1)*FF(:,1,h)'; 
SShj(h,:)=sum((Y-prey2(:,:,h)).^2); 
end 
SSh=(sum(SShj'))'; 
Qhj=ones(px-1,py)-PRESShj(2:px,:)./SShj(1:px-1,:); % 错位 
Qh=ones(px-1,1)-PRESSh(2:px,1)./SSh(1:px-1,1);