关于最大相关性最小冗余性(MRMR)的理论推导 http://t.csdnimg.cn/dzg2N

%%  主程序举例
all_f=all_feature;          %导入计算好特征向量矩阵all_feature
all_feature=mapminmax(all_f);           %统一化处理
label=[1*ones(1,30) 2*ones(1,30)]';     %贴标签
[mrmr,c]=mRMR(all_feature',label,20);   %调用mRMR子程序 
select_feature=all_f(c,:);                   

%%% 特征优选——最大相关最小冗余（MRMR）
%%%使用增量搜索算法来得到合适的特征集合
function [mrmr,slctFea] = mRMR(dataX, dataC, nSelect)
% 输入: 
%       dataX ，n-by-p：输入n*p特征向量矩阵，n个采样值，p个特征值；
%       dataC ，n-by-1：输入n*1向量，分类的结果（或贴好标签）；
%          nSelect  ：选择特征的个数，但是要小于或等于n个。

% 输出:
%       slctFea ，nSelect-by-1：输出n*1向量，优选的特征（排序）；
%                mrmr       : 相对应优选特征的值
%%  
%初始化
slctFea     = zeros(nSelect, 1);   %优选特征的零向量
dimF        = size(dataX, 2);      %特征的个数
remainFea   = (1:dimF);            %取前dimF个特征的序号

% 互信息（最大相关性），特征dataX与类别dataC之间，  维数为dim*1
Ixc = bundleMI(dataX, dataC);      %调用互信息子函数

% 选择第一个最大相关性的特征
[mrmr(1), idxF]          = max(Ixc);    %计算每一列中的最大元素mrmr(1)，以及他们所在的行索引idxF
slctFea(1)         = remainFea(idxF);   %存放第一个优选特征

% 删除已选定的特征
remainFea( idxF )  = [];
Ixc(       idxF )  = [];
dimRemain          = dimF-1;    % 剩余特征的个数

% 互信息（最小冗余性） ， 特征xi和特征xj之间
Ixx        = zeros(dimRemain, nSelect-1);
redunSum   = zeros(dimRemain, 1);

%% 使用增量搜索算法来得到合适的特征集合
for fea = 2 : nSelect    %去除第一个特征后
    
    Ixx(:, fea-1) = bundleMI( dataX(:,remainFea) , dataX(:,slctFea(fea-1)) );     %特征间互信息
    redunSum      = redunSum + Ixx(:,fea-1);

    % 选择具有增量mRMR的特性
    [mrmr(fea), idxF]   = max( Ixc - redunSum/(fea-1) );      %最大相关最小冗余特征集
    slctFea(fea)= remainFea(idxF);                            %%存放第fea个优选特征
   
    % 删除已选定的特征
    Ixx(idxF, : ) = [];           %特征与特征间的互信息
    redunSum(idxF) = [];
    remainFea(idxF) = [];
    Ixc( idxF) = [];              %特征与类别间的互信息
    dimRemain  = dimRemain - 1;   % 剩余特征的个数

end

end

%%% mRMR中的子程序，互信息
function bMI = bundleMI(X, y)
% 输入: 
%           X：维数为k的向量
%           y：维数为1的离散数据
%输出: 
%           bMI：维度为k*1,bMI(i) = I(X(:,i),y)
%%
k = size(X,2);       %向量X的维数
bMI = zeros(k,1);
for i = 1 : k
        bMI(i) = mutualInfoDis(X(:,i), y);    %调用另一子程序
end
end

%%% MI的子程序
function mi = mutualInfoDis(x, y)

n = length(x);      %向量x的长度
x = ceil(x(:));     %取大于或者等于指定表达式的最小整数
y = ceil(y(:));

%更新x和y，取值范围为1:max([x;y])
lo = min( [ x; y ] );
x  = x - lo + 1;
y  = y - lo + 1;
up = max( [ x; y ] );

%% 求熵
idx  = (1:n);
tabX = sparse(idx, x, 1, n, up, n);    %稀疏矩阵
                                       %由向量idx,x,1(三者维数一样)生成一个n*up的含有n个非零元素的稀疏矩阵tabX
                                       %idx 和 x 分别是矩阵中非零元素的行索引和列索引的向量。
                                       %1 是由对应的 (idx ,x) 对指定索引的非零值的向量。
                                       %最后一个n为：分配可以存储n个非零值的空间。                                     
tabY = sparse(idx, y, 1, n, up, n);
Pxy  = nonzeros(tabX'*tabY) / n;      %非零矩阵元素,概率p           
Px      = mean(tabX);
Py      = mean(tabY);
negHxy  = Pxy' * log(Pxy);             %x,y,联合熵   
negHx   = Px * log(Px)';             %x,熵，
negHy   = Py * log(Py)';             %y,熵，
%% 互信息
mi      = negHxy - negHx - negHy;     %互信息公式：H(x)+H(y)-H(x,y)
end