这是本人编写的PLS建模程序。
function [pls_P,pls_Q,pls_W,pls_B]=pls_model(X_M, Y_M, iter)
% reference: Paul Geladi and Bruce R. Kowalski. Partial Least-Square Regression: A Tutorial. Analytica Chimica Acta,1986, 185:19
% Fundation
% X=TP'+E=∑tp'+E
% Y=UQ'+F=∑uq'+F
% u=bt[pls_P,pls_Q,pls_W,pls_B]=pls_model(X_M, Y_M, iter)
% b=u't/t't
% Y=TBQ'+F
% input
% X_M N by K the independant block, must be processed with meaning and scaling
% Y_M N by P the dependant block, must be processed with meaning and scaling
% iter-- the number of PLS components extracted
% output
% pls_P the loadings for every pls component of X
% pls_Q the loadings for every pls componentof Y
% pls_W the weights
% pls_B the regression coefficients for the inner relation.
[N,K]=size(X_M);
[N,P]=size(Y_M);
if iter > min(N,K)
iter = min(N,K);
end
% start algrithm: PLS
for counter = 1:iter; % extract 3 PLS components
pls_u = Y_M(:,1); % get any column from Y
t_old_n = 1.0;% the ini of t old
t_new_n = 0.0;% the ini of t new
pls_t = 0;
while ( abs( t_old_n - t_new_n )>0.0001)
pls_t_old = pls_t;
% in the X block
pls_w = (pls_u' * X_M / ( pls_u' * pls_u ))';
pls_w = pls_w/norm(pls_w);
pls_t = X_M * pls_w;
% in the Y block
if(P>1)
pls_q = (pls_t' * Y_M /(pls_t'*pls_t))';
pls_q = pls_q/norm(pls_q);
pls_u = Y_M*pls_q;
else
pls_q=1;
end
% compare t,check for convergence
t_old_n = norm( pls_t_old );
t_new_n = norm( pls_t);
end
pls_p = (pls_t'*X_M/(pls_t'*pls_t))';
pls_p_n=norm(pls_p);
pls_p = pls_p/pls_p_n;
pls_t = pls_t*pls_p_n;
pls_w = pls_w*pls_p_n;
pls_b = pls_u'*pls_t/(pls_t'*pls_t);
X_M = X_M - pls_t*pls_p';
Y_M = Y_M - pls_b*pls_t*pls_q';
% pls_U(:,counter)=pls_u;
% pls_T(:,counter)=pls_t;
pls_P(:,counter)=pls_p;
pls_Q(:,counter)=pls_q;
pls_W(:,counter)=pls_w;
pls_B(:,counter)=pls_b;
end
% Now that correction is completed, start validation.
% Y_P = 0;
% for counter = 1:iter
% t_p = X_P*pls_W(:,counter);
% X_P=X_P-t_p*pls_P(:,counter)';
% Y_P =Y_P + pls_B(:,counter)*t_p*pls_Q(:,counter)';
% end
% for i = 1:m;
% Y_P(i,:)=Y_P(i,:).*std_Y+ave_Y_new;
% end
评论0