www.pudn.com > Successive-Projections-Algorithm.zip > spa.m, change:2007-03-30,size:4500b


function [var_sel,var_sel_phase2] = spa(Xcal,ycal,Xval,yval,m_min,m_max,autoscaling) 
 
% [var_sel,var_sel_phase2] = spa(Xcal,ycal,Xval,yval,m_min,m_max,autoscaling) --> Validation with a separate set 
% [var_sel,var_sel_phase2] = spa(Xcal,ycal,[],[],m_min,m_max,autoscaling) --> Cross-validation 
% 
% If m_min = [], the default value m_min = 1 is employed 
% If m_max = [], the default values m_max = min(N-1,K) (validation with a separate set) 
% or min(N-2,K) (cross-validation) are employed. 
% autoscaling --> 1 (yes) or 0 (no) 
 
autoscaling 
if ((autoscaling ~= 1) & (autoscaling ~= 0)) 
    error('Please choose whether or not to use autoscaling.') 
end 
 
N = size(Xcal,1); % Number of calibration objects 
K = size(Xcal,2); % Total number of variables 
 
if length(m_min) == 0, m_min = 1; end  
if length(m_max) == 0,  
    if size(Xval,1) > 0 
        m_max = min(N-1,K);  
    else 
        m_max = min(N-2,K);  
    end 
end  
 
if m_max > min(N-1,K) 
    error('m_max is too large !'); 
end 
 
% Phase 1: Projection operations for the selection of candidate subsets 
     
% The projections are applied to the columns of Xcal after 
% mean-centering and (optional) autoscaling 
     
if autoscaling == 1 
    normalization_factor = std(Xcal); 
else 
    normalization_factor = ones(1,K); 
end 
     
for k = 1:K 
    x = Xcal(:,k); 
    Xcaln(:,k) = (x - mean(x)) / normalization_factor(k); 
end 
 
SEL = zeros(m_max,K); 
 
h = waitbar(0,'Phase 1 (Projections). Please wait...'); 
loopStart = now; 
 
for k = 1:K 
    SEL(:,k) = projections_qr(Xcaln,k,m_max); 
    loopEnd = loopStart + (now-loopStart)*K/k; 
    waitbar(k/K,h,['Phase 1 ETC: ' datestr(loopEnd)]); 
end 
close(h); 
     
disp('Phase 1 (projections) completed !') 
     
% Phase 2: Evaluation of the candidate subsets according to the PRESS criterion 
 
PRESS = Inf*ones(m_max,K); 
h = waitbar(0,'Phase 2 (Evaluation of variable subsets). Please wait...'); 
warning off MATLAB:singularMatrix 
warning off MATLAB:nearlySingularMatrix 
loopStart = now; 
 
for k = 1:K 
    for m = m_min:m_max 
        var_sel = SEL(1:m,k); 
        [yhat,e] = validation(Xcal,ycal,Xval,yval,var_sel); 
        PRESS(m,k) = e'*e; 
    end 
    loopEnd = loopStart + (now-loopStart)*K/k; 
    waitbar(k/K,h,['Phase 2 ETC: ' datestr(loopEnd)]); 
end 
close(h); 
warning on MATLAB:singularMatrix 
warning on MATLAB:nearlySingularMatrix 
 
[PRESSmin,m_sel] = min(PRESS); 
[dummy,k_sel] = min(PRESSmin); 
 
var_sel_phase2 = SEL(1:m_sel(k_sel),k_sel); 
 
disp('Phase 2 (evaluation of variable subsets) completed !') 
 
% Phase 3: Final elimination of variables 
 
% Step 3.1: Calculation of the relevance index 
Xcal2 = [ones(N,1) Xcal(:,var_sel_phase2)];  
b = Xcal2\ycal; % MLR with intercept term 
std_deviation = std(Xcal2); 
relev = abs(b.*std_deviation'); 
relev = relev(2:end); % The intercept term is always included 
% Sorts the selected variables in decreasing order of "relevance" 
[dummy,index_increasing_relev] = sort(relev); % Increasing order 
index_decreasing_relev = index_increasing_relev(end:-1:1); % Decreasing order 
 
% Step 3.2: Calculation of PRESS values 
for i = 1:length(var_sel_phase2) 
    [yhat,e] = validation(Xcal,ycal,Xval,yval,var_sel_phase2(index_decreasing_relev(1:i)) ); 
    PRESS_scree(i) = e'*e; 
end 
RMSEP_scree = sqrt(PRESS_scree/length(e)); 
figure, grid, hold on 
plot(RMSEP_scree) 
xlabel('Number of variables included in the model'),ylabel('RMSE') 
 
% Step 3.3: F-test criterion 
PRESS_scree_min = min(PRESS_scree); 
alpha = 0.25; 
dof = length(e); % Number of degrees of freedom 
fcrit = finv(1-alpha,dof,dof); % Critical F-value 
PRESS_crit = PRESS_scree_min*fcrit; 
% Finds the minimum number of variables for which PRESS_scree 
% is not significantly larger than PRESS_scree_min  
i_crit = min(find(PRESS_scree < PRESS_crit));  
i_crit = max(m_min,i_crit); % The number of selected variables must be at least m_min 
 
var_sel = var_sel_phase2( index_decreasing_relev(1:i_crit) ); 
title(['Final number of selected variables: ' num2str(length(var_sel)) '  (RMSE = ' num2str(RMSEP_scree(i_crit)) ')']) 
 
% Indicates the selected point on the scree plot 
plot(i_crit,RMSEP_scree(i_crit),'s') 
 
disp('Phase 3 (final elimination of variables) completed !') 
 
% Presents the selected variables  
% in the first object of the calibration set 
figure,plot(Xcal(1,:));hold,grid 
plot(var_sel,Xcal(1,var_sel),'s') 
legend('First calibration object','Selected variables') 
xlabel('Variable index')