plstineng.m

MATLAB科学计算与工程分析源代码源程序3

X=[191 36 50;
189 37 52;
193 38 58;
162 35 62;
189 35 46;
182 36 56;
211 38 56;
167 34 60;
176 31 74;
154 33 56;
169 34 50;
166 33 52;
154 34 64;
247 46 50;
193 36 46;
202 37 62;
176 37 54;
157 32 52;
156 33 54;
138 33 68
];
Y=[5 162 60;
2 110 60;
12 101 101;
12 105 37;
13 155 58;
4 101 42;
8 101 38;
6 125 40;
15 200 40;
17 251 250;
17 120 38;
13 210 115;
14 215 105;
1 50 50 ;
6 70 31;
12 210 120;
4 60 25;
11 230 80;
15 225 73;
2 110 43];

[n,p]=size(X);
[n,q]=size(Y);

meanX=mean(X);%junzhi
varX=var(X);%fangcha
meanY=mean(Y);%junzhi
varY=var(Y);%fangcha


%%%%shu ju biao zhun hua
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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%omega1(:,1)=X0'*Y0/((Y0'*X0*X0'*Y0)^0.5)
[omega(:,1),t(:,1),pp(:,1),XX(:,:,1),rr(:,1),YY(:,:,1)]=tiquzhuyuan(X0,Y0);
%tt(:,1)=X0*omega1(:,1)
%ppp(:,1)=X0'*tt(:,1)/(tt(:,1)'*tt(:,1))
%XXX(:,:,1)=X0-tt(:,1)*ppp(:,1)'%%%di yi  ci zhu  yuan  

[omega(:,2),t(:,2),pp(:,2),XX(:,:,2),rr(:,2),YY(:,:,2)]=tiquzhuyuan(XX(:,:,1),YY(:,:,1));
%omega1(:,2)=XXX(:,:,1)'*Y0/((Y0'*XXX(:,:,1)*XXX(:,:,1)'*Y0)^0.5)
%tt(:,2)=XXX(:,:,1)*omega1(:,2)
%ppp(:,2)=XXX(:,:,1)'*tt(:,2)/(tt(:,2)'*tt(:,2))
%XXX(:,:,2)=XXX(:,:,1)-tt(:,2)*ppp(:,2)'  %%%di er ci  zhu yuan 

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,:);
  
    % for m=1:q
    %tttt=[ones(n-1,1),tt0];
   % [rrrrr,bint2,Residual2(:,m)]=regress(YY0(:,m),tttt);
    %rest(i,m)=Y0(i,m)-t(i,:)*rrrrr(2:3);
    %end
    
     expPRESS(i,:)=(Y0(i,:)-t(i,:)*inv((tt0'*tt0))*tt0'*YY0);
     for m=1:q
        PRESShj=PRESShj+expPRESS(i,m)^2;
      % PRESShj1=PRESShj+rest(i,m)^2;
   end
end
 sum1=sum(PRESShj);
 PRESSh=sum(sum1);
 
 for m=1:q
    %tttt=[ones(n,1),t(:,1)];
     %[r,bint2,Residual2(:,m)]=regress(Y0(:,m),tttt);
 
    for i=1:n
       % SShj(i,m)=(Residual2(i,m))^2;
       SShj(i,m)=YY(i,m,1)^2;
    end
 end
 sum2=sum(SShj);
 SSh=sum(sum2);

 Q=1-(PRESSh/SSh);

 k=3;
 %%%%%%%%%%%%%%%%  xun  huan,ti qu  zhu  yuan 
 while Q>0.0975
     [omega(:,k),t(:,k),pp(:,k),XX(:,:,k),rr(:,k),YY(:,:,k)]=tiquzhuyuan(XX(:,:,k-1),YY(:,:,k-1));
     %omega1(:,k)=XXX(:,:,k-1)'*Y0/((Y0'*XXX(:,:,k-1)*XXX(:,:,k-1)'*Y0)^0.5)
    %tt(:,k)=XXX(:,:,k-1)*omega1(:,k)
    %ppp(:,k)=XXX(:,:,k-1)'*tt(:,k)/(tt(:,k)'*tt(:,k))
    %XXX(:,:,k)=XXX(:,:,k-1)-tt(:,k)*ppp(:,k)' %%%di k ci  zhu yuan 
        
    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,:);
    
       %for m=1:q
   % tttt=[ones(n-1,1),tt00];
   % [rrrrr,bint2,Residual2(:,m)]=regress(YY0(:,m),tttt)
    %rest(i,m)=Y0(i,m)-t(i,:)*rrrrr(2:k+1);
   %end
           
     expPRESS(i,:)=(Y0(i,:)-t(i,:)*((tt00'*tt00)^(-1))*tt00'*YY0);
     for m=1:q
        PRESShj=PRESShj+expPRESS(i,m)^2;
        % PRESShj=PRESShj+rest(i,m)^2;
     end
 end

 
PRESShj;
sum1=sum(PRESShj);
 PRESSh=sum(sum1);
  
 for m=1:q
   % tttt=[ones(n,1),tt(:,1:(k-1))];
     %[r,bint2,Residual2(:,m)]=regress(Y0(:,m),tttt);
     for i=1:n
        % SShj(i,m)=(Residual2(i,m))^2;
        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%%%%%%%%% ti chu zhu yuan de ge shu 
%%%%%%%%%%%%%%         huan  yuan  xi  shu 

%omegaxing=ones(p,h);
%omegaxing(:,1)=r(2)*omega1(:,1);
%for i=2:(h)
  %for j=1:(i-1)
    %  omegaxingi =(eye(p)-omega1(:,j)*ppp(:,j)');
     % omegaxingii=eye(p);
     % omegaxingii=omegaxingii*omegaxingi;
  % end
     % omegaxing(:,i)=r(i+1)*omegaxingii*omega1(:,i);
%end

%beta1=sum(omegaxing,2);



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

%%%%%%%%  zui zhong de  xi  shu 
beta%%%%%%%%%%%%
%%%%%%%%%%%%

omega
omegaxing

%%%%%%%  ji  suan  xiang guan  xi  shu 
for i=1:h
    for j=1:q
        relation(i,j)=sum(prod(corrcoef(t(:,i),Y(:,j))))/2;
    end
end

relation

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Rd=relation.*relation;
RdYt=sum(Rd,2)/q
Rdtttt=sum(RdYt);
omega22=omega.*omega;

omega22*RdYt


VIP=((p/Rdtttt)*(omega22*RdYt)).^0.5;




clear