method_b2.m

Variable Reduction Testbench通过对变量进行相关性分析来实现减少变量的目的。

function [ranking, load_x, cut_off] = method_B2(X, is_corrmatrix, method)

% -------------------------------------------------------------------------
% this code is part of the 'Reduction Testbench' suite
% developed by A. Manganaro, R. Todeschini, A. Ballabio, D. Mauri
% 2006 - Milano Chemometrics and QSAR Research Group
% -------------------------------------------------------------------------
%
%
% [ranking, load_x, cut_off] = method_B2(X, is_corrmatrix, method)
%
% method_B2 uses Jolliffe's B2 algorithm to to evaluate an elimination ranking 
% of variables without any loss of data; it is suggested to delete
% until m variables are left, where m is calculated with the BS algorithm
% This routine also outputs a plot of the ranking
%
% Input:
% X = data set [n x p]  n objects, p variables
% is_corrmatrix = if set to 'y', X is seen as the dataset's correlation
%   matrix
% method = if set to 'i' enables the iterative version of B2
%
% Output:
% ranking = elimination ranking [1 x p] of the variables
% load_x = values [1 x p] of the EV% associated with the PC of each variable
% cut_off = number of variables suggested to be retained


echo off;

[n,p] = size(X);

if ( (p<2) | (n<2) )
    disp('Wrong matrix dimension - execution aborted');
    ranking=0; load_x=0; return;
end


% Sets the progressbar
screensize = [0 0 1 1];
width = screensize(3)/3;
height = screensize(4)/25;
left = screensize(3)/2 - width/2;
bottom = screensize(4)/2 - height/2;
progress_win = figure('Units','normalized','Position',[[left bottom] width height],...
                          'MenuBar','none','Resize','off','Name','Algorithm working...',...
                          'NumberTitle','off','WindowStyle','modal');
progress_axes = axes('Position',[0.02 0.15 0.96 0.70],'XLim',[0 1],'YLim',[0 1],'Box','on',...
                         'xtick',[],'ytick',[]);
progress_patch = patch('XData',[0 1 1 0],'YData',[0 0 1 1],'FaceColor',[1 0 0]);
drawnow;

p_tot = p;

set(progress_win,'Pointer','watch');


% Calculate the correlation matrix
if (is_corrmatrix)
    C = X;
else
    C = corrcoef(X);
end


% Calculate the number m of significative PC using the BS algorithm
[L,Lambda] = svd(C);
m = 0;
cur_ev = 1;
cur_dm = 0;
while (m<p) & (cur_ev>cur_dm)
    m = m + 1;
    cur_ev = (Lambda(m,m) / trace(Lambda)) * 100;
    tmpsum=0;
    for idx = m:p
        tmpsum = tmpsum + 1/idx;
    end
    cur_dm = (100/p)*tmpsum;
end
m = m - 1;


if (method=='i')

    % Iterative method
    
    rank_idx = 1;
    for idx=1:p
        var_id(idx)=idx;
    end


    while (p>1)
    
        % Updates the progressbar
        set(progress_win,'Name',['Algorithm working... ',int2str(p),' variables left']);
        set(progress_patch,'XData',[0 p/p_tot p/p_tot 0],'Facecolor',...
            [1 (1-p/p_tot) 0]);
        drawnow;

        % Eigenvalues and eigenvectors of C
        [L,Lambda] = svd(C);
    
        % Finds the most representative variable of the last PC
        L(:,p) = abs(L(:,p));
        max_var = find(L(:,p)==max(L(:,p)));

        % Updates ranking and load_x
        ranking(rank_idx) = var_id(max_var(1));
        load_x(rank_idx) = (Lambda(p,p)/sum(diag(Lambda)))*100;
    
        % Deletes from the dataset the variable selected
        if (max_var(1)==1)
            var_id = var_id(:,2:end);
            C = C(2:end,2:end);
        elseif (max_var(1)==p)
            var_id = var_id(:,1:(p-1));
            C = C(1:(p-1),1:(p-1));
        else
            var_id = [var_id(:,1:max_var(1)-1),var_id(:,max_var(1)+1:p)];
            C = C([1:(max_var(1)-1) (max_var(1)+1):p],[1:(max_var(1)-1) (max_var(1)+1):p]);
        end
    
        % Updates indexes
        rank_idx = rank_idx + 1;
        p = p-1;
    
    end

    % Adds to ranking the last variable
    ranking = [ranking, var_id];
    load_x = [load_x,   (Lambda(1,1)/sum(diag(Lambda)))*100];

else
    
    % Non-iterative method

    ranking = [];
    rank_idx = 1;
    
    while (p>0)
    
        % Updates the progressbar
        set(progress_win,'Name',['Algorithm working... ',int2str(p),' variables left']);
        set(progress_patch,'XData',[0 p/p_tot p/p_tot 0],'Facecolor',...
            [1 (1-p/p_tot) 0]);
        drawnow;

        cur_pc(:,1) = [1:length(L(:,p))]';
        cur_pc(:,2) = abs(L(:,p));
        
        cur_pc = sortrows(cur_pc,2);

        max_var = length(cur_pc);
        
        while (max_var>1)
            flag = 0;
            for i=1:length(ranking)
                if (ranking(i)==cur_pc(max_var,1))
                    flag=1;
                    break;
                end
            end
            
            if (flag==0)
                break;
            end
            
            max_var = max_var - 1;
        end

        % Updates ranking and selected
        ranking(rank_idx) = cur_pc(max_var,1);
        load_x(rank_idx) = (Lambda(p,p)/sum(diag(Lambda)))*100;
        
        % Updates indexes
        rank_idx = rank_idx + 1;
        p = p-1;
    
    end
  
end
    
% Closes the progressbar
set(progress_win,'Pointer','arrow');
close(progress_win);

cut_off = m;


%%%% Fig_1: graph of resulting ranking %%%%

fig_1 = figure('Name','Variables elimination ranking - B2 method');
title('Variables ranking by B2 method');
xlabel('Variables ranking');
ylabel('EV% of the correspondent PC');
ylim([0 max(load_x)+(10*max(load_x)/100)]);
xlim([0 length(load_x)+1]);

hold on;

line([(length(load_x)-m+0.6) (length(load_x)-m+0.6)],ylim,'LineStyle','--','Color','r');
plot(load_x,'o-b','MarkerFaceColor','b');
for i=1:length(ranking)
    text(i,load_x(i)++(5*max(load_x)/100),[' ',num2str(ranking(i))]);
end
text((length(load_x)-m+0.7),(5*max(load_x)/100),['variables = ',num2str(m)]);

hold off;


echo on;