\max _{x_j\in{X-S_{m-1}}}[I(x_j;c)-\frac{1}{m-1}\sum_{x_i\in S_{m-1}}I(x_i;x_j)]

该算法计算复杂度为O(|S| M)

3. 算法实现

function [fea] = mrmr_miq_d(d, f, K)

% function [fea] = mrmr_miq_d(d, f, K)

% d-输入N*M矩阵，N个采样值，M个特征，特征向量

% f-输入的N*1向量，分类结果

% K-选择的葛铮的个数

% MIQ scheme according to MRMR

%

% By Hanchuan Peng

% April 16, 2003

%

bdisp=0;

nd = size(d,2);%特征个数

nc = size(d,1);%采样个数

t1=cputime;

for i=1:nd

   t(i) = mutualinfo(d(:,i), f);%分别计算每个特征与分类f之间的互信息

end

fprintf('calculate the marginal dmi costs %5.1fs.\n', cputime-t1);

[tmp, idxs] = sort(-t);%按照计算的互信息降序排列，tmp存储的是降序排列的值，idxs为原来的序号

fea_base = idxs(1:K);%取前k个特征的序号，赋值给fea_base

fea(1) = idxs(1);%将互信息最大的特征赋值给fea(1)

KMAX = min(1000,nd); %如果特征多于1000个，则首先选择前1000个

idxleft = idxs(2:KMAX);%获得剩余特征的序号

k=1;

if bdisp==1

fprintf('k=1 cost_time=(N/A) cur_fea=%d #left_cand=%d\n', ...

      fea(k), length(idxleft));

end

for k=2:K

   t1=cputime;

   ncand = length(idxleft);%获得剩余特征的个数

   curlastfea = length(fea);%获得当前已选择的特征数

   for i=1:ncand

      t_mi(i) = mutualinfo(d(:,idxleft(i)), f); %计算剩余特征与分类f之间的互信息，相当于公式中的D

      mi_array(idxleft(i),curlastfea) = getmultimi(d(:,fea(curlastfea)), d(:,idxleft(i)));%计算每个特征之间的互信息

      c_mi(i) = mean(mi_array(idxleft(i), :)); %求特征之间互信息的平均值，相当于公式中的R

   end

%   [tmp, fea(k)] = max(t_mi(1:ncand) ./ c_mi(1:ncand));

   [tmp, fea(k)] = max(t_mi(1:ncand) ./ (c_mi(1:ncand) + 0.01));

    %此处与公式中不一样，此处用的是D/R的最大值作为标准，+0.01是为了避免c_mi为0的情况

   tmpidx = fea(k); fea(k) = idxleft(tmpidx); idxleft(tmpidx) = [];

    %选择第k个特征

   if bdisp==1

   fprintf('k=%d cost_time=%5.4f cur_fea=%d #left_cand=%d\n', ...

      k, cputime-t1, fea(k), length(idxleft));

   end

end

return;

%===================================== 

function c = getmultimi(da, dt) 

for i=1:size(da,2)

   c(i) = mutualinfo(da(:,i), dt);

end