my_pls_pred3.m

基于交叉检验的PLS算法

function [OptVecNum,OptTheta,SSEPred,SSE]=my_Pls_Pred3(X,Y)%采用pls建立预测模型，输出测试数据中的实测输出值和模型预测值


Xbuild=X(1:400,1:12);%建模训练输入数据
Xtest=X(401:600,1:12);%测试输入数据
Ybuild=Y(1:400,1);%建模训练输出数据
Ytest=Y(401:600,1);%测试输出数据

%调用函数unify，输入输出数据归一化
meanX=mean(Xbuild);
stdX=std(Xbuild);
meanY=mean(Ybuild);
stdY=std(Ybuild);
Xbuild=unify(Xbuild,meanX,stdX);
Xtest=unify(Xtest,meanX,stdX);
Ybuild=unify(Ybuild,meanY,stdY);
Ytest=unify(Ytest,meanY,stdY);

%将建模训练数据分为四组
Xbuild1=Xbuild(1:100,:);
Xbuild2=Xbuild(101:200,:);
Xbuild3=Xbuild(201:300,:);
Xbuild4=Xbuild(301:400,:);
Ybuild1=Ybuild(1:100,1);
Ybuild2=Ybuild(101:200,1);
Ybuild3=Ybuild(201:300,1);
Ybuild4=Ybuild(301:400,1);

%利用交叉检验法确定特征向量个数
SSE=[];
Theta=[];
for VecNum=1:1:12
    
    theta1=pls([Xbuild1;Xbuild2;Xbuild3],[Ybuild1;Ybuild2;Ybuild3],VecNum);%调用pls函数
    SSE1=norm((Ybuild4-Xbuild4*theta1),2);%预测误差平方和SSE=sum of squared errors
    
    theta2=pls([Xbuild1;Xbuild2;Xbuild4],[Ybuild1;Ybuild2;Ybuild4],VecNum);%调用pls函数
    SSE2=norm((Ybuild3-Xbuild3*theta2),2);%预测误差平方和SSE=sum of squared errors
    
    theta3=pls([Xbuild1;Xbuild4;Xbuild3],[Ybuild1;Ybuild4;Ybuild3],VecNum);%调用pls函数
    SSE3=norm((Ybuild2-Xbuild2*theta3),2);%预测误差平方和SSE=sum of squared errors
    
    theta4=pls([Xbuild4;Xbuild2;Xbuild3],[Ybuild4;Ybuild2;Ybuild3],VecNum);%调用pls函数
    SSE4=norm((Ybuild1-Xbuild1*theta4),2);%预测误差平方和SSE=sum of squared errors
    
    Theta=[Theta (theta1+theta2+theta3+theta4)/4];
    SSE(VecNum)=(SSE1+SSE2+SSE3+SSE4)/4;
end
%取最优特征向量个数
OptVecNum=1;
for i=2:1:12
    if SSE(i)<SSE(OptVecNum)
        OptVecNum=i;
    end
end  
%取最优模型参数thet
OptTheta=Theta(:,OptVecNum);
%做SSE图
%subplot(1,2,1)
%plot(SSE);
%grid on;
%xlabel('number of principal component');ylabel('SSE');
%做实测输出和预测输出图
%subplot(1,2,2)
Ypred=Xtest*OptTheta;
SSEPred=norm((Ytest-Ypred),2);%检验SSE
plot(Ytest,'r-');
hold on;
plot(Ypred,'b:');
title('Actual(red -) and Predicted(blue --) Output');
xlabel('Sample Number');ylabel('Output');


function [X]=unify(X,meanX,stdX)%数据归一化

[n,m]=size(X);
for j=1:m
    for i=1:n
        X(i,j)=(X(i,j)-meanX(1,j))/stdX(1,j);
    end
end


function [theta]=pls(X,Y,VecNum)%partial least squares,VecNum＝number of eigenvector

%初始化
E=X;F=Y;
T=[];P=[];Q=[];U=[];W=[];
u=Y(:,1);
[n,m]=size(X);
t=ones(n,1);
%外部变换
for h=1:VecNum
    tempt=zeros(n,1);
    while norm((t-tempt),2)>1e-6
        tempt=t;
        tempw=u'*E/(u'*u);
        w=tempw';
        w=w/norm(w,2);
        t=E*w/(w'*w);
        tempc=t'*F/(t'*t);
        c=tempc';
        c=c/norm(c,2);    
        u=F*c/(c'*c);
    end
    Optc=fminsearch(@(optc) myfun(optc,t,u,n),[0;inv(t'*t)*u'*t;0]);%利用fminsearch函数寻优
    uNew=[ones(n,1) t t.^2]*Optc;
    tempwNew=uNew'*E/(uNew'*uNew);
    wNew=tempwNew';
    wNew=wNew/norm(wNew,2);
    tNew=E*wNew/(wNew'*wNew);    
    t=tNew;
    u=uNew;
    w=wNew;
    p=E'*tNew/(tNew'*tNew);
    p=p/norm(p,2);
    q=F'*uNew/(uNew'*uNew);
    q=q/norm(q,2);
    E=E-tNew*p';
    F=F-uNew*q';
    T=[T tNew];P=[P p];Q=[Q q];U=[U uNew];W=[W wNew];
end

theta=W*inv(P'*W)*inv(T'*T)*T'*Y;