基于互信息的特征选择算法MATLAB实现

uqHZau

/ a+ s) y0 R% Z2 c在概率论和信息论中，两个随机变量的互信息（Mutual Information，简称MI）或转移信息（transinformation）是变量间相互依赖性的量度。不同于相关系数，互信息并不局限于实值随机变量，它更加一般且决定着联合分布 p(X,Y) 和分解的边缘分布的乘积 p(X)p(Y) 的相似程度。互信息(Mutual Information)是度量两个事件集合之间的相关性(mutual dependence)。互信息最常用的单位是bit。

7 M- `2 M% O& J' U. Y. X4 w互信息的定义
2 c2 B$ I- ?! t正式地，两个离散随机变量 X 和 Y 的互信息可以定义为：
9 I. W9 a( g( O; ^
* X# N6 G! j( w  I, e4 D% O4 ?* p其中 p(x,y) 是 X 和 Y 的联合概率分布函数，而p(x)和p(y)分别是 X 和 Y 的边缘概率分布函数。
& f/ l5 g" u9 i) c/ N
% e8 `' m9 `- W* y
6 j6 x2 [) e/ J# K1 z
# M7 e3 K- ?8 {5 j8 L! M/ ]在连续随机变量的情形下，求和被替换成了二重定积分：


5 W% h6 [0 b7 S2 b& B* a
$ r& W9 b5 e2 `. b0 L- K5 V: W! f其中 p(x,y) 当前是 X 和 Y 的联合概率密度函数，而p(x)和p(y)分别是 X 和 Y 的边缘概率密度函数。
) W: C2 k" d  J& J# c& N
( I; _& t  A1 W' v  Y7 S互信息量I(xi;yj)在联合概率空间P(XY)中的统计平均值。 平均互信息I(X;Y)克服了互信息量I(xi;yj)的随机性,成为一个确定的量。如果对数以 2 为基底，互信息的单位是bit。
. J" F& ?( w& \5 {6 y' I7 c
0 I) X; Q( ?  d$ L  j) j直观上，互信息度量 X 和 Y 共享的信息：它度量知道这两个变量其中一个，对另一个不确定度减少的程度。例如，如果 X 和 Y 相互独立，则知道 X 不对 Y 提供任何信息，反之亦然，所以它们的互信息为零。在另一个极端，如果 X 是 Y 的一个确定性函数，且 Y 也是 X 的一个确定性函数，那么传递的所有信息被 X 和 Y 共享：知道 X 决定 Y 的值，反之亦然。因此，在此情形互信息与 Y（或 X）单独包含的不确定度相同，称作 Y（或 X）的熵。而且，这个互信息与 X 的熵和 Y 的熵相同。（这种情形的一个非常特殊的情况是当 X 和 Y 为相同随机变量时。）
, ^, ]) m' |0 s0 F
+ b, R  Z# F! V, J2 r: U互信息是 X 和 Y 联合分布相对于假定 X 和 Y 独立情况下的联合分布之间的内在依赖性。于是互信息以下面方式度量依赖性：I(X; Y) = 0 当且仅当 X 和 Y 为独立随机变量。从一个方向很容易看出：当 X 和 Y 独立时，p(x,y) = p(x) p(y)，因此：

7 U5 ]( |1 z$ f, G5 l
" c5 g, R5 [; |
此外，互信息是非负的（即 I(X;Y) ≥ 0; 见下文），而且是对称的（即 I(X;Y) = I(Y;X)）。

5 c: N$ A3 p- O$ u$ T3 h3 B7 @互信息特征选择算法的步骤
3 j( O  L# z* n" N

• ①划分数据集
• ②利用互信息对特征进行排序
• ③选择前n个特征利用SVM进行训练
• ④在测试集上评价特征子集计算错误率
  - }3 d  X% o% q) i5 N

缺点
此种特征选择方法是最大化特征与分类变量之间的相关度，就是选择与分类变量拥有最高相关度的前k个变量。但是，在特征选择中，单个好的特征的组合并不能增加分类器的性能，因为有可能特征之间是高度相关的，这就导致了特征变量的冗余。

代码
注意使用的数据集是dlbcl，大概五千多维，可以从UCI上下载，最终选择前100特征进行训练。
( D9 ^% f' O! N4 U4 n  b
% O* I5 P3 p4 r, u" L主函数代码：
! w5 D$ ?5 e/ k$ y+ h; h9 o1 q

• clear all
• close all
• clc;
• [X_train,Y_train,X_test,Y_test] = divide_dlbcl();
• Y_train(Y_train==0)=-1;
• Y_test(Y_test==0)=-1;
• % number of features
• numF = size(X_train,2);
• [ ranking , w] = mutInfFS( X_train, Y_train, numF );
• k = 100; % select the Top 2 features
• svmStruct = svmtrain(X_train(:,ranking(1:k)),Y_train,'showplot',true);
• C = svmclassify(svmStruct,X_test(:,ranking(1:k)),'showplot',true);
• err_rate = sum(Y_test~= C)/size(X_test,1); % mis-classification rate
• conMat = confusionmat(Y_test,C); % the confusion matrix
• fprintf('\nAccuracy: %.2f%%, Error-Rate: %.2f \n',100*(1-err_rate),err_rate);
  

   
' a" P8 g1 F+ r
mutInfFS.m

6 W7 M. b, E& s' }/ p# S1 e: @

• function [ rank , w] = mutInfFS( X,Y,numF )
• rank = [];
• for i = 1:size(X,2)
•     rank = [rank; -muteinf(X(:,i),Y) i];
• end;
• rank = sortrows(rank,1);
• w = rank(1:numF, 1);
• rank = rank(1:numF, 2);
• end
  + o( i8 Y; z! ^( S1 d

% P0 n$ x$ u+ x, `; q
. b- a9 v$ ~' u: M2 o/ Rmuteinf.m

' N2 {6 v, f! N) u' B, p4 n3 t

• function info = muteinf(A, Y)
• n = size(A,1);%实例数量
• Z = [A Y];%所有实例的维度值及标签
• if(n/10 > 20)
•     nbins = 20;
• else
•     nbins = max(floor(n/10),10);%设置区间的个数
• end;
• pA = hist(A, nbins);%min(A)到max(A)划分出nbins个区间出来，求每个区间的概率
• pA = pA ./ n;%除以实例数量
• i = find(pA == 0);
• pA(i) = 0.00001;%不能使某一区间的概率为0
• od = size(Y,2);%一个维度
• cl = od;
• %下面是求实例不同标签的的概率值，也就是频率
• if(od == 1)
•     pY = [length(find(Y==+1)) length(find(Y==-1))] / n;
•     cl = 2;
• else
•     pY = zeros(1,od);
•     for i=1：od
•         pY(i) = length(find(Y==+1));
•     end;
•     pY = pY / n;
• end;
• p = zeros(cl,nbins);
• rx = abs(max(A) - min(A)) / nbins;%每个区间长度
• for i = 1:cl
•     xl = min(A);%变量的下界
•     for j = 1:nbins
•         if(i == 2) && (od == 1)
•             interval = (xl <= Z(:,1)) & (Z(:,2) == -1);
•         else
•             interval = (xl <= Z(:,1)) & (Z(:,i+1) == +1);
•         end;
•         if(j < nbins)
•             interval = interval & (Z(:,1) < xl + rx);
•         end;
•         %find(interval)
•         p(i,j) = length(find(interval));
•         if p(i,j) == 0 % hack!
•             p(i,j) = 0.00001;
•         end
•         xl = xl + rx;
•     end;
• end;
• HA = -sum(pA .* log(pA));%计算当前维度的信息熵
• HY = -sum(pY .* log(pY));%计算标签的信息熵
• pA = repmat(pA,cl,1);
• pY = repmat(pY',1,nbins);
• p = p ./ n;
• info = sum(sum(p .* log(p ./ (pA .* pY))));
• info = 2 * info ./ (HA + HY);%计算互信息
  

   
" n" B' Q  [. u) G3 y
前100个特征的效果：
& B/ h3 N& e" o
4 d2 ]; z& z# X3 P  kAccuracy: 86.36%, Error-Rate: 0.14
1 v5 l0 v! m# \6 B
# J$ m. \/ {" w) a选择前两个特征进行训练（压缩率接近100%,把上述代码中的K设为2即可）的二维图：
/ \1 v4 X+ e* c6 @6 r, O1 z! v

, X' y/ R. r$ j
Accuracy: 75.00%, Error-Rate: 0.25

ExxNEN

基于互信息的特征选择算法MATLAB实现