所谓偏最小二乘法

1、 所谓偏最小二乘法，就是指在做基于最小二乘法的线性回归分析之前，对数据集进行主成分分析降维，下面的源码是没有删减的，GreenSim团队免费提供您使用，转载请注明GreenSim团队（/greensim）。function y5,e1,e2=PLS(X,Y,x,y,p,q)% 偏最小二乘回归的通用程序% 注释以“基于近红外光谱分析的汽油组分建模”为例，但本程序的适用范围绝不仅限于此% GreenSim团队原创作品（/greensim）% 输入参数列表% X 校正集光谱矩阵，nk的矩阵，n个样本，k个波长%

2、Y 校正集浓度矩阵，nm的矩阵，n个样本，m个组分% x 验证集光谱矩阵% y 验证集浓度矩阵% p X的主成分的个数，最佳取值需由其它方法确定% q Y的主成分的个数，最佳取值需由其它方法确定% 输出参数列表% y5 x对应的预测值（y为真实值）% e1 预测绝对误差，定义为e1=y5-y% e2 预测相对误差，定义为e2=|(y5-y)/y|% 第一步：对X,x,Y,y进行归一化处理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

3、(j)=max(Xx(:,j); Xx(:,j)=(Xx(:,j)-xmin(j)/(xmax(j)-xmin(j);endymin=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);endX1=Xx(1:n,:);x1=Xx(n+1):end,:);Y1=Yy(1:n,:);y1=Yy(n+1):end,:); 第二步：分别提取X1和Y1的p和q个主成分，并将X1,x1,Y1,y1映射到主成分空间CX,

4、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; % 第三步：对X2和Y2进行线性回归B=regress(Y2,X2,0.05);%第三个输入参数是显著水平，可以调整 % 第四步：将x2带入模型得到预测值y3y3=x2*B; % 第五步：将y3进行“反主成分变换”得到y4y4=y3*pinv(CY);% 第六步：将y4反归一化得到y5for j=1:m y5(:,j)=(ymax(j)-ymin(j)*y4(:,j)+ymin(j);en

5、d% 第七步：计算误差e1=y5-y;e2=abs(y5-y)./y);function MD,ERROR,PRESS,SECV,SEC=ExtraSim1(X,Y)% 基于PLS方法的进一步仿真分析% 功能一：计算MD值，以便于发现奇异样本% 功能二：计算各种p取值情况下的ERROR,PRESS,SECV,SEC值，以确定最佳输入变量个数%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=si

6、ze(XX,1);for p=1:pmax disp(p); Err1=zeros(1,N);%绝对误差 Err2=zeros(1,N);%相对误差 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,:

7、); 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偏最小二乘法 matlab程序函数：function T,P,U,Q,B,W = pl

8、s (X,Y,tol2)% PLS Partial Least Squares Regrassion% T,P,U,Q,B,Q = pls(X,Y,tol) performs particial least squares regrassion% between the independent variables, X and dependent Y as% X = T*P + E;% Y = U*Q + F = T*B*Q + F1;% Inputs:% X data matrix of independent variables 自变量% Y data matrix of dependen

9、t variables因变量% tol the tolerant of convergence (defaut 1e-10) 收敛度% % Outputs: Y = X*B + E, X = T*P + E0, and Y = U*Q + F0, U = T*B + F1% T score matrix of X 得分矩阵% P loading matrix of X 载荷矩阵% U score matrix of Y% Q loading matrix of Y% B matrix of regression coefficient 矩阵回归系数% W weight matrix of X

10、权重矩阵% Using the PLS model, for new X1, Y1 can be predicted as% Y1 = (X1*P)*B*Q = X1*(P*B*Q)% or% Y1 = X1*(W*inv(P*W)*inv(T*T)*T*Y)% Without Y provided, the function will return the principal components as% X = T*P + E% Example: taken from Geladi, P. and Kowalski, B.R., An example of 2-block% predict

11、ive partial least-squares regression with simulated data,% Analytica Chemica Acta, 185(1996) 19-32.%x=4 9 6 7 7 8 3 2;6 15 10 15 17 22 9 4;8 21 14 23 27 36 15 6;10 21 14 13 11 10 3 4; 12 27 18 21 21 24 9 6; 14 33 22 29 31 38 15 8;16 33 22 19 15 12 3 6; 18 39 26 27 25 26 9 8;20 45 30 35 35 40 15 10;y

12、=1 1;3 1;5 1;1 3;3 3;5 3;1 5;3 5;5 5;% leave the last sample for testN=size(x,1);x1=x(1:N-1,:); %取x中18行的所有数据y1=y(1:N-1,:);x2=x(N,:); %取x中的第9行数据y2=y(N,:);% normalization 归一化xmean=mean(x1);xstd=std(x1);ymean=mean(y1);ystd=std(y);X=(x1-xmean(ones(N-1,1),:)./xstd(ones(N-1,1),:);Y=(y1-ymean(ones(N-1,1),:

13、)./ystd(ones(N-1,1),:);% PLS modelT,P,U,Q,B,W=pls(X,Y);% Prediction and erroryp = (x2-xmean)./xstd * (P*B*Q);fprintf(Prediction error: %g/n,norm(yp-(y2-ymean)./ystd);% By Yi Cao at Cranfield University on 2nd Febuary 2008% Reference:% Geladi, P and Kowalski, B.R., Partial Least-Squares Regression: A

14、% Tutorial, Analytica Chimica Acta, 185 (1986) 1-7.% Input checkerror(nargchk(1,3,nargin);error(nargoutchk(0,6,nargout);if nargin2 Y=X;endtol = 1e-10;if nargintol2 & k tol w = X*u; w = w/norm(w); t = t1; t1 = X*w; q = Y*t1; q = q/norm(q); u = Y*q; end % update p based on t t=t1; p=X*t/(t*t); pnorm=n

15、orm(p); p=p/pnorm; t=t*pnorm; w=w*pnorm; % regression and residuals b = u*t/(t*t); X = X - t*p; Y = Y - b*t*q; % save iteration results to outputs: k=k+1; T(:,k)=t; P(:,k)=p; U(:,k)=u; Q(:,k)=q; W(:,k)=w; B(k,k)=b; % uncomment the following line if you wish to see the convergence% disp(norm(Y)endT(:

16、,k+1:end)=;P(:,k+1:end)=;U(:,k+1:end)=;Q(:,k+1:end)=;W(:,k+1:end)=;B=B(1:k,1:k); 例子：% Principal Component Analysis and Partial Least Squares% Principal Component Analysis (PCA) and Partial Least Squares (PLS) are% widely used tools. This code is to show their relationship through the% Nonlinear Iter

17、ative PArtial Least Squares (NIPALS) algorithm.% The Eigenvalue and Power Method% The NIPALS algorithm can be derived from the Power method to solve the% eigenvalue problem. Let x be the eigenvector of a square matrix, A,% corresponding to the eignvalue s:% $Ax=sx$% Modifying both sides by A iterati

18、vely leads to% $Akx=skx$% Now, consider another vectro y, which can be represented as a linear% combination of all eigenvectors:% $y=/sum_in b_ix_i=Xb$% where% $X=/leftx_1/,/,/, /cdots/,/,/, x_n /right$% and% $b = /leftb_1/,/,/, /cdots/,/,/, b_n /rightT$% Modifying y by A gives% $Ay=AXb=XSb$% Where

19、S is a diagnal matrix consisting all eigenvalues. Therefore, for% a large enough k,% $Aky=XSkb/approx /alpha x_1$% That is the iteration will converge to the direction of x_1, which is the% eigenvector corresponding to the eigenvalue with the maximum module.% This leads to the following Power method

20、 to solve the eigenvalue problem.A=randn(10,5);% sysmetric matrix to ensure real eigenvaluesB=A*A;%find the column which has the maximum normdum,idx=max(sum(A.*A);x=A(:,idx);%storage to judge convergencex0=x-x;%convergence toleranttol=1e-6;%iteration if not convergedwhile norm(x-x0)tol %iteration to

21、 approach the eigenvector direction y=A*x; %normalize the vector y=y/norm(y); %save previous x x0=x; %x is a product of eigenvalue and eigenvector x=A*y;end% the largest eigen value corresponding eigenvector is ys=x*x;% compare it with those obtained with eigV,D=eig(B);d,idx=max(diag(D);v=V(:,idx);d

22、isp(d-s)% v and y may be different in signsdisp(min(norm(v-y),norm(v+y)% The NIPALS Algorithm for PCA% The PCA is a dimension reduction technique, which is based on the% following decomposition:% $X=TPT+E$% Where X is the data matrix (m x n) to be analysed, T is the so called% score matrix (m x a),

23、P the loading matrix (n x a) and E the residual.% For a given tolerance of residual, the number of principal components, a,% can be much smaller than the orginal variable dimension, n. % The above power algorithm can be extended to get T and P by iteratively% subtracting A (in this case, X) by x*y (

24、in this case, t*p) until the% given tolerance satisfied. This is the so called NIPALS algorithm.% The data matrix with normalizationA=randn(10,5);meanx=mean(A);stdx=std(A);X=(A-meanx(ones(10,1),:)./stdx(ones(10,1),:);B=X*X;% allocate T and PT=zeros(10,5);P=zeros(5);% tol for convergencetol=1e-6;% to

25、l for PC of 95 percenttol2=(1-0.95)*5*(10-1);for k=1:5 %find the column which has the maximum norm dum,idx=max(sum(X.*X); t=A(:,idx); %storage to judge convergence t0=t-t; %iteration if not converged while norm(t-t0)tol %iteration to approach the eigenvector direction p=X*t; %normalize the vector p=

26、p/norm(p); %save previous t t0=t; %t is a product of eigenvalue and eigenvector t=X*p; end %subtracing PC identified X=X-t*p; T(:,k)=t; P(:,k)=p; if norm(X)tol2 break endendT(:,k+1:5)=;P(:,k+1:5)=;S=diag(T*T);% compare it with those obtained with eigV,D=eig(B);D,idx=sort(diag(D),descend);D=D(1:k);V=

27、V(:,idx(1:k);fprintf(The number of PC: %i/n,k);fprintf(norm of score difference between EIG and NIPALS: %g/n,norm(D-S);fprintf(norm of loading difference between EIG and NIPALS: %g/n,norm(abs(V)-abs(P);% The NIPALS Algorithm for PLS% For PLS, we will have two sets of data: the independent X and depe

28、ndent% Y. The NIPALS algorithm can be used to decomposes both X and Y so that% $X=TPT+E,/,/,/,/,Y=UQT+F,/,/,/,/,U=TB$% The regression, U=TB is solved through least sequares whilst the% decompsition may not include all components. That is why the approach is% called partial least squares. This algori

29、thm is implemented in the PLS% function.% Example: Discriminant PLS using the NIPALS Algorithm% From Chiang, Y.Q., Zhuang, Y.M and Yang, J.Y, Optimal Fisher% discriminant analysis using the rank decomposition, Pattern Recognition,% 25 (1992), 101-111.% Three classes data, each has 50 samples and 4 v

30、ariables.x1=5.1 3.5 1.4 0.2; 4.9 3.0 1.4 0.2; 4.7 3.2 1.3 0.2; 4.6 3.1 1.5 0.2;. 5.0 3.6 1.4 0.2; 5.4 3.9 1.7 0.4; 4.6 3.4 1.4 0.3; 5.0 3.4 1.5 0.2; . 4.4 2.9 1.4 0.2; 4.9 3.1 1.5 0.1; 5.4 3.7 1.5 0.2; 4.8 3.4 1.6 0.2; . 4.8 3.0 1.4 0.1; 4.3 3.0 1.1 0.1; 5.8 4.0 1.2 0.2; 5.7 4.4 1.5 0.4; . 5.4 3.9 1

31、.3 0.4; 5.1 3.5 1.4 0.3; 5.7 3.8 1.7 0.3; 5.1 3.8 1.5 0.3; . 5.4 3.4 1.7 0.2; 5.1 3.7 1.5 0.4; 4.6 3.6 1.0 0.2; 5.1 3.3 1.7 0.5; . 4.8 3.4 1.9 0.2; 5.0 3.0 1.6 0.2; 5.0 3.4 1.6 0.4; 5.2 3.5 1.5 0.2; . 5.2 3.4 1.4 0.2; 4.7 3.2 1.6 0.2; 4.8 3.1 1.6 0.2; 5.4 3.4 1.5 0.4; . 5.2 4.1 1.5 0.1; 5.5 4.2 1.4

32、0.2; 4.9 3.1 1.5 0.2; 5.0 3.2 1.2 0.2; . 5.5 3.5 1.3 0.2; 4.9 3.6 1.4 0.1; 4.4 3.0 1.3 0.2; 5.1 3.4 1.5 0.2; . 5.0 3.5 1.3 0.3; 4.5 2.3 1.3 0.3; 4.4 3.2 1.3 0.2; 5.0 3.5 1.6 0.6; . 5.1 3.8 1.9 0.4; 4.8 3.0 1.4 0.3; 5.1 3.8 1.6 0.2; 4.6 3.2 1.4 0.2; . 5.3 3.7 1.5 0.2; 5.0 3.3 1.4 0.2;x2=7.0 3.2 4.7 1

33、.4; 6.4 3.2 4.5 1.5; 6.9 3.1 4.9 1.5; 5.5 2.3 4.0 1.3; . 6.5 2.8 4.6 1.5; 5.7 2.8 4.5 1.3; 6.3 3.3 4.7 1.6; 4.9 2.4 3.3 1.0; . 6.6 2.9 4.6 1.3; 5.2 2.7 3.9 1.4; 5.0 2.0 3.5 1.0; 5.9 3.0 4.2 1.5; . 6.0 2.2 4.0 1.0; 6.1 2.9 4.7 1.4; 5.6 2.9 3.9 1.3; 6.7 3.1 4.4 1.4; . 5.6 3.0 4.5 1.5; 5.8 2.7 4.1 1.0;

34、 6.2 2.2 4.5 1.5; 5.6 2.5 3.9 1.1; . 5.9 3.2 4.8 1.8; 6.1 2.8 4.0 1.3; 6.3 2.5 4.9 1.5; 6.1 2.8 4.7 1.2; . 6.4 2.9 4.3 1.3; 6.6 3.0 4.4 1.4; 6.8 2.8 4.8 1.4; 6.7 3.0 5.0 1.7; . 6.0 2.9 4.5 1.5; 5.7 2.6 3.5 1.0; 5.5 2.4 3.8 1.1; 5.5 2.4 3.7 1.0; . 5.8 2.7 3.9 1.2; 6.0 2.7 5.1 1.6; 5.4 3.0 4.5 1.5; 6.

35、0 3.4 4.5 1.6; . 6.7 3.1 4.7 1.5; 6.3 2.3 4.4 1.3; 5.6 3.0 4.1 1.3; 5.5 2.5 5.0 1.3; . 5.5 2.6 4.4 1.2; 6.1 3.0 4.6 1.4; 5.8 2.6 4.0 1.2; 5.0 2.3 3.3 1.0; . 5.6 2.7 4.2 1.3; 5.7 3.0 4.2 1.2; 5.7 2.9 4.2 1.3; 6.2 2.9 4.3 1.3; . 5.1 2.5 3.0 1.1; 5.7 2.8 4.1 1.3;x3=6.3 3.3 6.0 2.5; 5.8 2.7 5.1 1.9; 7.1

36、 3.0 5.9 2.1; 6.3 2.9 5.6 1.8; . 6.5 3.0 5.8 2.2; 7.6 3.0 6.6 2.1; 4.9 2.5 4.5 1.7; 7.3 2.9 6.3 1.8; . 6.7 2.5 5.8 1.8; 7.2 3.6 6.1 2.5; 6.5 3.2 5.1 2.0; 6.4 2.7 5.3 1.9; . 6.8 3.0 5.5 2.1; 5.7 2.5 5.0 2.0; 5.8 2.8 5.1 2.4; 6.4 3.2 5.3 2.3; . 6.5 3.0 5.5 1.8; 7.7 3.8 6.7 2.2; 7.7 2.6 6.9 2.3; 6.0 2.

37、2 5.0 1.5; . 6.9 3.2 5.7 2.3; 5.6 2.8 4.9 2.0; 7.7 2.8 6.7 2.0; 6.3 2.7 4.9 1.8; . 6.7 3.3 5.7 2.1; 7.2 3.2 6.0 1.8; 6.2 2.8 4.8 1.8; 6.1 3.0 4.9 1.8; . 6.4 2.8 5.6 2.1; 7.2 3.0 5.8 1.6; 7.4 2.8 6.1 1.9; 7.9 3.8 6.4 2.0; . 6.4 2.8 5.6 2.2; 6.3 2.8 5.1 1.5; 6.1 2.6 5.6 1.4; 7.7 3.0 6.1 2.3; . 6.3 3.4 5.6 2.4; 6.4 3.1 5.5 1.8; 6.0 3.0 4.8 1.8; 6.9 3.1 5.4 2.1; . 6.7 3.1 5.6 2.4; 6.9 3.1 5.1 2.3; 5.8 2.7 5.1 1.9; 6.8 3.2 5.9 2.3; . 6.7 3.3 5.7 2.5; 6.7 3.0 5.2 2.3; 6.3 2.5 5.0 1.9; 6.5 3.0 5.2 2.0; . 6.2 3.4 5.4 2.3; 5.9 3.0 5.1 1.8;%Split data