www.pudn.com > spa.rar > SSPPAA.m, change:2016-04-20,size:2817b


clc 
clear 
close all 
load data_chaye2.mat 
Xcal=newcy(1:30,1:405); 
ycal=newcy(1:30,5799); 
Xval=newcy(31:50,1:405); 
yval=newcy(31:50,5799); 
m_min=1; 
m_max=200; 
 
N = size(Xcal,1);  
K = size(Xcal,2);  
 
 
% Phase 1: Projection operations for the selection of candidate subsets 
     
normalization_factor = std(Xcal); 
     
for k = 1:K 
    x = Xcal(:,k); 
    Xcaln(:,k) = (x - mean(x)) / normalization_factor(k); 
end 
 
 
SEL = zeros(m_max,K); 
for k = 1:K 
    SEL(:,k) = projection(Xcaln,k,m_max); 
end 
 
     
% Phase 2: Evaluation of the candidate subsets according to the PRESS criterion 
 
PRESS = Inf*ones(m_max,K); 
 
for k = 1:K 
    for m = m_min:m_max 
        var_sel = SEL(1:m,k); 
        [yhat1,e1] = validation(Xcal,ycal,Xval,yval,var_sel); 
        PRESS(m,k) = e1'*e1; 
    end 
end 
 
 
[PRESSmin,m_sel] = min(PRESS);  %找出每列最小值,返回所在行(波段数) 
[dummny,k_sel] = min(PRESSmin); %找到矩阵里最小值 
 
%第k_sel波段为初始波段时最佳,波段数目为m_sel(k_sel) 
var_sel_phase2 = SEL(1:m_sel(k_sel),k_sel);  
 
 
% 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'); 
relev1 = relev(2:end); 
% Sorts the selected variables in decreasing order of "relevance" 
[dummy,index_increasing_relev] = sort(relev1); 
index_decreasing_relev = index_increasing_relev(end:-1:1); 
% Step 3.2: Calculation of PRESS values 
for i = 1:length(var_sel_phase2) 
    [yhat2,e2] = validation(Xcal,ycal,Xval,yval,var_sel_phase2(index_decreasing_relev(1:i)) ); 
    PRESS_scree(i) = e2'*e2; 
end 
RMSEP_scree = sqrt(PRESS_scree/length(e2)); 
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(e2); % 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  
zx=find(PRESS_scree < PRESS_crit); 
i_crit = min(zx);  
i_crit = max(m_min,i_crit);  
 
var_sel1 = 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') 
 
 
% Presents the selected variables  
% in the first object of the calibration set 
figure,plot(Xcal(6,:));hold,grid 
plot(var_sel1,Xcal(6,var_sel1),'s') 
legend('First calibration object','Selected variables') 
xlabel('Variable index') 
save boduan_shang var_sel1