【原創(chuàng)】偏最小二乘法的Matlab源碼
所謂偏最小二乘法，就是指在做基于最小二乘法的線性回歸分析之前，對(duì)數(shù)據(jù)集進(jìn)行主成分分析降維，下面的源碼是沒有刪減的，GreenSim團(tuán)隊(duì)免費(fèi)提供您使用，轉(zhuǎn)載請(qǐng)注明GreenSim團(tuán)隊(duì)（http://blog.sina.com.cn/greensim）。

function [y5,e1,e2]=PLS(X,Y,x,y,p,q)
%% 偏最小二乘回歸的通用程序
%  注釋以“基于近紅外光譜分析的汽油組分建模”為例，但本程序的適用范圍絕不僅限于此
%  GreenSim團(tuán)隊(duì)原創(chuàng)作品（http://blog.sina.com.cn/greensim）
%% 輸入?yún)?shù)列表
% X        校正集光譜矩陣，n×k的矩陣，n個(gè)樣本，k個(gè)波長
% Y        校正集濃度矩陣，n×m的矩陣，n個(gè)樣本，m個(gè)組分
% x        驗(yàn)證集光譜矩陣
% y        驗(yàn)證集濃度矩陣
% p        X的主成分的個(gè)數(shù)，最佳取值需由其它方法確定
% q        Y的主成分的個(gè)數(shù)，最佳取值需由其它方法確定
%% 輸出參數(shù)列表
% y5       x對(duì)應(yīng)的預(yù)測(cè)值（y為真實(shí)值）
% e1       預(yù)測(cè)絕對(duì)誤差，定義為e1=y5-y
% e2       預(yù)測(cè)相對(duì)誤差，定義為e2=|(y5-y)/y|

%% 第一步：對(duì)X,x,Y,y進(jìn)行歸一化處理
[n,k]=size(X);
m=size(Y,2);
Xx=[X;x];
Yy=[Y;y];
xmin=zeros(1,k);
xmax=zeros(1,k);
for j=1:k
    xmin(j)=min(Xx(:,j));
    xmax(j)=max(Xx(:,j));
    Xx(:,j)=(Xx(:,j)-xmin(j))/(xmax(j)-xmin(j));
end
ymin=zeros(1,m);
ymax=zeros(1,m);
for j=1:m
    ymin(j)=min(Yy(:,j));
    ymax(j)=max(Yy(:,j));
    Yy(:,j)=(Yy(:,j)-ymin(j))/(ymax(j)-ymin(j));
end
X1=Xx(1:n,:);
x1=Xx((n+1):end,:);
Y1=Yy(1:n,:);
y1=Yy((n+1):end,:);

%% 第二步：分別提取X1和Y1的p和q個(gè)主成分，并將X1,x1,Y1,y1映射到主成分空間
[CX,SX,LX]=princomp(X1);
[CY,SY,LY]=princomp(Y1);
CX=CX(:,1:p);
CY=CY(:,1:q);
X2=X1*CX;
Y2=Y1*CY;
x2=x1*CX;
y2=y1*CY;

%% 第三步：對(duì)X2和Y2進(jìn)行線性回歸
B=regress(Y2,X2,0.05);%第三個(gè)輸入?yún)?shù)是顯著水平，可以調(diào)整

%% 第四步：將x2帶入模型得到預(yù)測(cè)值y3
y3=x2*B;

%% 第五步：將y3進(jìn)行“反主成分變換”得到y(tǒng)4
y4=y3*pinv(CY);

%% 第六步：將y4反歸一化得到y(tǒng)5
for j=1:m
    y5(:,j)=(ymax(j)-ymin(j))*y4(:,j)+ymin(j);
end

%% 第七步：計(jì)算誤差
e1=y5-y;
e2=abs((y5-y)./y);

function [MD,ERROR,PRESS,SECV,SEC]=ExtraSim1(X,Y)
%% 基于PLS方法的進(jìn)一步仿真分析
%% 功能一：計(jì)算MD值，以便于發(fā)現(xiàn)奇異樣本
%% 功能二：計(jì)算各種p取值情況下的ERROR,PRESS,SECV,SEC值，以確定最佳輸入變量個(gè)數(shù)
%  GreenSim團(tuán)隊(duì)原創(chuàng)作品（http://blog.sina.com.cn/greensim）
%%
[n,k]=size(X);
m=size(Y,2);
pmax=n-1;
q=m;
ERROR=zeros(1,pmax);
PRESS=zeros(1,pmax);
SECV=zeros(1,pmax);
SEC=zeros(1,pmax);
XX=X;
YY=Y;
N=size(XX,1);
for p=1:pmax
    disp(p);
    Err1=zeros(1,N);%絕對(duì)誤差
    Err2=zeros(1,N);%相對(duì)誤差
    for i=1:N
        disp(i);
        if i==1
            x=XX(1,:);
            y=YY(1,:);
            X=XX(2:N,:);
            Y=YY(2:N,:);
        elseif i==N
            x=XX(N,:);
            y=YY(N,:);
            X=XX(1:(N-1),:);
            Y=YY(1:(N-1),:);
        else
            x=XX(i,:);
            y=YY(i,:);
            X=[XX(1:(i-1),:);XX((i+1):N,:)];
            Y=[YY(1:(i-1),:);YY((i+1):N,:)];
        end
        [y5,e1,e2]=PLS(X,Y,x,y,p,q);
        Err1(i)=e1;
        Err2(i)=e2;
    end
    ERROR(p)=sum(Err2)/N;
    PRESS(p)=sum(Err1.^2);
    SECV(p)=sqrt(PRESS(p)/n);
    SEC(p)=sqrt(PRESS(p)/(n-p));
end
%%
[CX,SX,LX]=princomp(X);
S=SX(:,1:p);
MD=zeros(1,n);
for j=1:n
    s=S(j,:);
    MD(j)=(s')*(inv(S'*S))*(s);
end