No science, No life.

もっと科学を身近に

相互情報量を計算してみた

「脳の情報表現」という本にある「相互情報量 (mutual information)」を計算してみました。

脳の情報表現―ニューロン・ネットワーク・数理モデル

脳の情報表現―ニューロン・ネットワーク・数理モデル

  • 作者:
  • 出版社/メーカー: 朝倉書店
  • 発売日: 2002/03
  • メディア: 単行本

エントロピーは、ある確率変数 X のもつ不確かさを表す量で、下のように定義されます。

 H(X) = \sum -P(x) \log P(x)

2変数の同時確率の場合は、

 H(X1, X2) = \sum -P(x1, x2) \log P(x1, x2)

です。そして、相互情報量は、

 I(X1, X2) = H(X1) + H(X2) - H(X1, X2)

と定義されます。この相互情報量という値は、一方の確率変数 X1 の値を知ったときに、もうひとつの確率変数 X2 についての不確かさの減少分を指しています。

この本の20ページ目にある下の例について、実際に計算してみました。

X1 と X2 が相関している場合。相互情報量 I(X1,X2) は 1.08 ビットとなりました (片方の X1 の値を知ると、1.08 ビット分だけ、X2 の値の不確かさが減る)。
f:id:neocortex:20190616233827p:plain

X1, X2 の分布は上と同じですが、片方をシャッフルして独立にした場合。相互情報量は低下して、0.10 ビットとなりました (片方の X1 の値を知っても、もう一方の X2 についての情報はほとんど得られない)。
f:id:neocortex:20190616233842p:plain

Matlab コードは下の通りです。

function mi()

xaxis = -3.9:0.2:3.9;
nbin = numel(xaxis);
bin = xaxis(2)-xaxis(1);

% 2d gaussian random variables
X = mvnrnd([0 0],[0.8 0.7; 0.7 0.8],5000);

% shuffling
if 0
    idx = randperm(size(X,1));
    X = [X(:,1) X(idx,2)];
end

% probability distribution function
pdf1 = hist(X(:,1),xaxis); pdf1 = pdf1/sum(pdf1); 
pdf2 = hist(X(:,2),xaxis); pdf2 = pdf2/sum(pdf2);
pdf12 = zeros(nbin,nbin);
for ii=1:nbin
    for jj=1:nbin
        pdf12(ii,jj) = sum(xaxis(ii)-bin/2<X(:,1) & X(:,1)<=xaxis(ii)+bin/2 &...
            xaxis(jj)-bin/2<X(:,2) & X(:,2)<=xaxis(jj)+bin/2);
    end
end
pdf12 = pdf12./sum(sum(pdf12));

% entoropy
H1 = getEntropy(pdf1);
H2 = getEntropy(pdf2);
H12= getEntropy(pdf12);

% mutual information
I12 = H1+H2-H12;

% plot
figure(1); clf
subplot(3,3,[4 8])
scatter(X(:,1),X(:,2),'.')
axis([-4 4 -4 4])
title(sprintf('H(x1,x2) = %.2f, I(x1,x2) = %.2f',H12,I12))
xlabel('x1');
ylabel('x2')

subplot(3,3,[1 2])
bar(xaxis,pdf1,1)
xlim([-4 4])
title(sprintf('H(x1) = %.2f',H1))
ylabel('P(x1)')

subplot(3,3,[6 9])
barh(xaxis,pdf2,1)
ylim([-4 4])
title(sprintf('H(x2) = %.2f',H2))
xlabel('P(x2)')

end

function H = getEntropy(pdf)

% linearlize
pdf = pdf(:);

% normalize (just in case)
pdf = pdf./sum(pdf);

% remove 0
pdf(pdf==0) = [];

% entropy 
H = sum(-pdf.*log2(pdf));

end