[projects] example code for mutual information

projects/project_mutualinfo/code/mi.m

@@ -0,0 +1,8 @@
+function I = mi(nxy)
+    pxy = nxy / sum(nxy(:));
+    px = sum(nxy, 2) / sum(nxy(:));
+    py = sum(nxy, 1) / sum(nxy(:));
+    pi = pxy .* log2(pxy./(px*py));
+    pi(nxy == 0) = 0.0;
+    I = sum(pi(:));
+end

projects/project_mutualinfo/code/mutualinfo.m

@@ -0,0 +1,90 @@
+%% load data:
+x = load('../data/decisions.mat');
+presented = x.presented;
+reported = x.reported;
+
+%% plot data:
+figure()
+plot(presented, 'ob', 'markersize', 10, 'markerfacecolor', 'b');
+hold on;
+plot(reported, 'or', 'markersize', 5, 'markerfacecolor', 'r');
+hold off
+ylim([0.5, 2.5])
+
+p1 = sum(presented == 1);
+p2 = sum(presented == 2);
+r1 = sum(reported == 1);
+r2 = sum(reported == 2);
+figure()
+bar([p1, p2, r1, r2]);
+set(gca, 'XTickLabel', {'p1', 'p2', 'r1', 'r2'});
+
+%% histogram:
+nxy = zeros(2, 2);
+for x = [1, 2]
+    for y = [1, 2]
+        nxy(x, y) = sum((presented == x) & (reported == y));
+    end
+end
+figure()
+bar3(nxy)
+set(gca, 'XTickLabel', {'p1', 'p2'});
+set(gca, 'YTickLabel', {'r1', 'r2'});
+
+%% normalized histogram:
+pxy = nxy / sum(nxy(:));
+figure()
+imagesc(pxy)
+
+px = sum(nxy, 2) / sum(nxy(:));
+py = sum(nxy, 1) / sum(nxy(:));
+
+%% mutual information:
+miv = mi(nxy);
+
+%% permutation:
+np = 10000;
+mis = zeros(np, 1);
+for k = 1:np
+    ppre = presented(randperm(length(presented)));
+    prep = reported(randperm(length(reported)));
+    pnxy = zeros(2, 2);
+    for x = [1, 2]
+        for y = [1, 2]
+            pnxy(x, y) = sum((ppre == x) & (prep == y));
+        end
+    end
+    mis(k) = mi(pnxy);
+end
+alpha = sum(mis>miv)/length(mis);
+fprintf('signifikance: %g\n', alpha);
+bins = [0.0:0.025:0.4];
+hist(mis, bins)
+hold on;
+plot([miv, miv], [0, np/10], '-r')
+hold off;
+xlabel('MI')
+ylabel('Count')
+
+%% maximum MI:
+n = 100000;
+pxs = [0:0.01:1.0];
+mis = zeros(length(pxs), 1);
+for k = 1:length(pxs)
+    p = rand(n, 1);
+    nxy = zeros(2, 2);
+    nxy(1, 1) = sum(p<pxs(k));
+    nxy(2, 2) = length(p) - nxy(1, 1);
+    mis(k) = mi(nxy);
+%nxy(1, 2) = 0;
+%nxy(2, 1) = 0;
+%mi(nxy)
+end
+figure();
+plot(pxs, mis);
+hold on;
+plot([px(1), px(1)], [0, 1], '-r')
+hold off;
+xlabel('p(x=1)')
+ylabel('Max MI=Entropy')
+