relieff.m

Feature Selection using matlab

function [ResultMat,FeaturesWeightsOrdered,FeaturesIndexOrdered,...
  OptimumFeatureSet] = ReliefF(DatasetToUse,ErrorEstMethod,handles)

%=============== Load The Patterns ===============================

% NPatterns: The number of Patterns
% KFeatures: The number of features
% CClasses : The number of features

% Patterns, features and Targets in a single matrix of 
% NPatterns X (KFeatures + 1) dimensionality. 
% The additional feature column is the Targets. 
% Patterns: FLOAT numbers in [0,1]
% Targets: INTEGER in {1,2,...,C}, where C the number of classes. 

global NPatterns KFeatures Patterns Targets CClasses Prior 
global StopByUser SetOfClasses

%================= Data Load ======================================
[Patterns, Targets] = DataLoadAndPreprocess(DatasetToUse);
 
[NPatterns, KFeatures] = size(Patterns);
CClasses = max(Targets);

for IndexClass =1:CClasses
    Prior(IndexClass) = sum(Targets == IndexClass)/ NPatterns;
end

%============ Settings of ReliefF =================================

NCorePatterns = 250; % floor(NPatterns/5); % User defined 1...NPatterns
NHits         = 10;                  % USer defined 
PercTest      = 10;                 % For CCR
NRepThres     = 50;
GammaParam    = 0.025;
ConfMatSwitch = 0;
SetOfClasses  = 1:CClasses;
NDc = NPatterns*(1-1/PercTest)/CClasses;

%============ Initialization ======================================
FeatureWeightsTab = zeros(KFeatures,1);
CurrentFeatureSet = [];
ResultMat         = [];
OptimumFeatureSet = [];
OptimumLowLimMahal = 0;
CCR                  = zeros(NCorePatterns,1);  
InfoLoss             = zeros(NCorePatterns,1);
LowLimitMahalInfLoss = zeros(NCorePatterns,1);
%============ Begin ReliefF========================================
for IndexCorePattern = 1:NCorePatterns
    CorePatternIndex = round(NPatterns*rand); % Select rand a pat
    
    if CorePatternIndex == 0
        CorePatternIndex = 1;
    end
    IndexCoreClass   = Targets(CorePatternIndex);
    CorePattern      = Patterns(CorePatternIndex,:);
    SetNoCoreClass   = SetOfClasses;
    SetNoCoreClass(IndexCoreClass) = [];
    IndexMisses=FindNearMisses(CorePattern, IndexCoreClass, NHits);
    IndexHits = FindNearHits(CorePattern, IndexCoreClass, NHits);
    for IndexFeature = 1:KFeatures
        FeatureWeightsTab(IndexFeature) = ...
                  FeatureWeightsTab(IndexFeature) ...
                       -1/NCorePatterns/NHits*...
                          sum(abs(CorePattern(IndexFeature)-...
                                Patterns(IndexHits,IndexFeature)));
                            
        for IndexNoCoreClass = SetNoCoreClass
           FeatureWeightsTab(IndexFeature) = ...
             FeatureWeightsTab(IndexFeature) + ...
              Prior(IndexNoCoreClass)/(1-Prior(IndexCoreClass))/...
                   NCorePatterns/NHits*...
                   sum(abs(CorePattern(IndexFeature) - ...
                     Patterns(IndexMisses(:,IndexNoCoreClass),...
                                                   IndexFeature)));
        end
    end
    [FeaturesWeightsOrdered, FeaturesIndexOrdered] = ...
                              sort( FeatureWeightsTab, 'descend');

    IndexImportFeats  = find(FeaturesWeightsOrdered > 0);
    if length(IndexImportFeats) > NDc - 1                          
    IndexImportFeats  = find(FeaturesWeightsOrdered > ...
                                  0.7*max(FeaturesWeightsOrdered));
    end
    
    CurrentFeatureSet = FeaturesIndexOrdered(IndexImportFeats);
    if length(CurrentFeatureSet) > NDc - 1
        CurrentFeatureSet = CurrentFeatureSet(1:NDc-1);
    end
    Dim               = length(CurrentFeatureSet);
    
    [CCR(IndexCorePattern), ConfMat, lowCL, upCL] = ...
            BayesClassMVGaussPDFs(Patterns(:,CurrentFeatureSet),...
                   Targets, PercTest, ErrorEstMethod, NRepThres,...
                                     GammaParam, 0, ConfMatSwitch);

    InfoLoss(IndexCorePattern) = CalcInfoLoss(Dim, floor(NDc),...
                                                   ErrorEstMethod);
       
    LowLimitMahalInfLoss(IndexCorePattern) = ...
           CCR(IndexCorePattern) - InfoLoss(IndexCorePattern)...
                              *(CCR(IndexCorePattern)-1/CClasses); 
    
    ResultMat(end+1,1:4) = [IndexCorePattern ...
     CCR(IndexCorePattern) LowLimitMahalInfLoss(IndexCorePattern) ...
     length(IndexImportFeats)];

    ResultMat
     
     
     
    % If best store it                               
    if LowLimitMahalInfLoss(IndexCorePattern)>OptimumLowLimMahal 
        OptimumFeatureSet = CurrentFeatureSet;
        OptimumStepIndx   = IndexCorePattern;
        OptimumLowLimMahal=LowLimitMahalInfLoss(IndexCorePattern);
    end
                                 
    %----------------------- Plot Module --------------------------
    if ~isempty(handles)
        StrList = ...
          num2str([FeaturesIndexOrdered FeaturesWeightsOrdered]);
        set(findobj(gcf,'Tag','ListSelFeats'), 'String', StrList);
        
        axes(handles.YelLinesAxes);
        axis([0 NPatterns 0 KFeatures]); axis manual
        hold on
        CountLines = 0;
        HYelLines = zeros(1,length(CurrentFeatureSet));
        for IndexFeature = CurrentFeatureSet'
            CountLines = CountLines + 1;
            if (NPatterns > KFeatures)
                HYelLines(CountLines) = plot([0 NPatterns+2],...
                                       IndexFeature*ones(1,2),'y');
            else
                HYelLines(CountLines)  =...
                 plot(IndexFeature*ones(1,2),[0 NPatterns+2], 'y');
            end
        end
        drawnow        
  
        axes(handles.FeatSelCurve);
        plot(1:IndexCorePattern, CCR(1:IndexCorePattern), 'b');
        hold on
        plot(1:IndexCorePattern, ...
                 LowLimitMahalInfLoss(1:IndexCorePattern), 'r');
        xlabel('# of core patterns');
        drawnow
        delete(HYelLines);
    else
        subplot(2,1,1)
        plot( ResultMat(1:IndexCorePattern,2),'b.-');
        hold on
        plot( ResultMat(1:IndexCorePattern,3), 'r.-');
        drawnow
        title([DatasetToUse ' ' ErrorEstMethod ' ReliefF']);
        subplot(2,1,2)
        hist(FeaturesWeightsOrdered,100)
%        axis([-0.02 0.02 0 100])
    end
    %------------------- End plot modulo --------------------------    
    if StopByUser
        StopByUser = 0;
        return
    end
     %-------------------------------------------------------------
end  % End NCorePatterns Repetition  


if ~isempty(handles)
    axes(handles.YelLinesAxes);
    for IndexFeature = OptimumFeatureSet'
        CountLines = CountLines + 1;
        hold on
        if (NPatterns > KFeatures)
            HYelLines(CountLines) = plot([0 NPatterns+2],...
                                       IndexFeature*ones(1,2),'y');
        else
            HYelLines(CountLines)  =...
                plot(IndexFeature*ones(1,2),[0 NPatterns+2], 'y');
        end
    end
    set(findobj(gcf,'Tag','ListSelFeats'), 'String', ...
                                                OptimumFeatureSet);
axes(handles.FeatSelCurve);
end

subplot(2,1,1)
hold on
plot(OptimumStepIndx*ones(1,2), [0 OptimumLowLimMahal], 'r');
StopByUser = 0;
return
%================ Calculate Negative Term =========================
function IndexHits = FindNearHits(CorePattern,IndexCoreClass,NHits)

global Patterns Targets 
IndexCoreClassPatterns = find(Targets == IndexCoreClass);
NIndexCoreClassPatterns = length(IndexCoreClassPatterns);

[A,B] = sort(sum(abs(repmat(CorePattern,...
          NIndexCoreClassPatterns,1) ... 
                      - Patterns(IndexCoreClassPatterns,:)),2));
IndexHits = B(2:NHits+1);
%================== Calculate Positive Term =======================
function IndexMisses =FindNearMisses(CorePattern,IndexCoreClass,...
                                                             NHits)

global Patterns Targets SetOfClasses CClasses
IndexMisses    = zeros(NHits,CClasses);
SetNoCoreClass = SetOfClasses;
SetNoCoreClass(IndexCoreClass) = [];
for IndexNoCoreClass = SetNoCoreClass
    IndexNoCoreClassPatterns = find(Targets == IndexNoCoreClass);
    NIndexNoCoreClassPatterns = length(IndexNoCoreClassPatterns);
    [A,B] = sort(sum(abs(...
        repmat(CorePattern,NIndexNoCoreClassPatterns,1)-...
                    Patterns(IndexNoCoreClassPatterns,:)),2));
    IndexMisses(:,IndexNoCoreClass) = B(1:NHits);
end
%==================================================================