%% 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')