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

clc
clear
close all
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