1. 程式人生 > >基於互資訊的特徵選擇演算法MATLAB實現

基於互資訊的特徵選擇演算法MATLAB實現

在概率論和資訊理論中,兩個隨機變數的互資訊(Mutual Information,簡稱MI)或轉移資訊(transinformation)是變數間相互依賴性的量度。不同於相關係數,互資訊並不侷限於實值隨機變數,它更加一般且決定著聯合分佈 p(X,Y) 和分解的邊緣分佈的乘積 p(X)p(Y) 的相似程度。互資訊(Mutual Information)是度量兩個事件集合之間的相關性(mutual dependence)。互資訊最常用的單位是bit。

互資訊的定義
正式地,兩個離散隨機變數 X 和 Y 的互資訊可以定義為:

其中 p(x,y) 是 X 和 Y 的聯合概率分佈函式,而p(x)和p(y)分別是 X 和 Y 的邊緣概率分佈函式。
這裡寫圖片描述


在連續隨機變數的情形下,求和被替換成了二重定積分:
這裡寫圖片描述
其中 p(x,y) 當前是 X 和 Y 的聯合概率密度函式,而p(x)和p(y)分別是 X 和 Y 的邊緣概率密度函式。

互資訊量I(xi;yj)在聯合概率空間P(XY)中的統計平均值。 平均互資訊I(X;Y)克服了互資訊量I(xi;yj)的隨機性,成為一個確定的量。如果對數以 2 為基底,互資訊的單位是bit。

直觀上,互資訊度量 X 和 Y 共享的資訊:它度量知道這兩個變數其中一個,對另一個不確定度減少的程度。例如,如果 X 和 Y 相互獨立,則知道 X 不對 Y 提供任何資訊,反之亦然,所以它們的互資訊為零。在另一個極端,如果 X 是 Y 的一個確定性函式,且 Y 也是 X 的一個確定性函式,那麼傳遞的所有資訊被 X 和 Y 共享:知道 X 決定 Y 的值,反之亦然。因此,在此情形互資訊與 Y(或 X)單獨包含的不確定度相同,稱作 Y(或 X)的熵。而且,這個互資訊與 X 的熵和 Y 的熵相同。(這種情形的一個非常特殊的情況是當 X 和 Y 為相同隨機變數時。)

互資訊是 X 和 Y 聯合分佈相對於假定 X 和 Y 獨立情況下的聯合分佈之間的內在依賴性。於是互資訊以下面方式度量依賴性:I(X; Y) = 0 當且僅當 X 和 Y 為獨立隨機變數。從一個方向很容易看出:當 X 和 Y 獨立時,p(x,y) = p(x) p(y),因此:
這裡寫圖片描述
此外,互資訊是非負的(即 I(X;Y) ≥ 0; 見下文),而且是對稱的(即 I(X;Y) = I(Y;X))。

更多互資訊內容請訪問:http://www.omegaxyz.com/2018/08/02/mi/

互資訊特徵選擇演算法的步驟
①劃分資料集
②利用互資訊對特徵進行排序
③選擇前n個特徵利用SVM進行訓練
④在測試集上評價特徵子集計算錯誤率
程式碼


注意使用的資料集是dlbcl,大概五千多維,可以從UCI上下載,最終選擇前100特徵進行訓練。

主函式程式碼:

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);

mutInfFS.m

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

muteinf.m

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);%計算互資訊

前100個特徵的效果:

Accuracy: 86.36%, Error-Rate: 0.14

選擇前兩個特徵進行訓練(壓縮率接近100%,把上述程式碼中的K設為2即可)的二維圖:
這裡寫圖片描述
Accuracy: 75.00%, Error-Rate: 0.25

更多內容訪問omegaxyz.com
網站所有程式碼採用Apache 2.0授權
網站文章採用知識共享許可協議BY-NC-SA4.0授權
© 2018 • OmegaXYZ-版權所有 轉載請註明出處