% enee324StatsBootEx

K = 6;
SX = 1:K;
rolls = [6, 2, 3, 2, 3, 6, 6, 3, 1, 4, 1, 6, 5, 2, 5, 1, 6, ...
         4, 2, 4, 6, 4, 5, 6, 5, 5, 5, 4, 2, 1, 6, 3, 6, 4, ...
         6, 6, 6, 3, 2, 2, 3, 4, 6, 6, 1, 4, 4, 5, 6, 1, 5, ...
         6, 6, 6, 6, 4, 6, 1, 4, 4, 2, 2, 2, 4, 3, 5, 2, 4, ...
         2, 5, 6, 4, 6, 2, 3, 6, 2, 4, 2, 5, 3, 6, 5, 2, 5, ...
         1, 1, 2, 6, 6, 2, 6, 3, 1, 1, 5, 2, 4, 1, 1 ];
n = length(rolls)
rollsByValue = hist(rolls,K)
m = 1/n*sum(rolls)
s2 = 1/(n-1)*sum((rolls-m).^2)
s = sqrt(s2)
se2 = s2/n
se = sqrt(se2)

contrast = max(rollsByValue)/min(rollsByValue)
s2PnP = 1/n*sum((rolls-m).^2)

B = 5000

rollsByValueB = nan(B,K);
mB = nan(B,1);
contrastB = nan(B,1);

for b =1:B
    indexTmp = ceil(n*rand(n,1));
    rollsTmp = rolls(indexTmp);
    rollsByValueB(b,:) = hist(rollsTmp,K);
    mB(b) = 1/n*sum(rollsTmp);
    contrastB(b) = max(rollsByValueB(b,:))/min(rollsByValueB(b,:));
end

figure(1),hist(mB,20)
mBoot = mean(mB)
seBoot = sqrt(var(mB))
mBootBias = mBoot - m % quite small, so ignorable
line([m m],ylim, 'Color','r');
line([mBoot mBoot],ylim, 'Color','m','LineWidth',2); % basically same as m
line([mBoot mBoot]-se,ylim, 'Color','k','LineWidth',2);
line([mBoot mBoot]+se,ylim, 'Color','k','LineWidth',2);
line([mBoot mBoot]-seBoot,ylim, 'Color','c','LineWidth',2); % basically same as se
line([mBoot mBoot]+seBoot,ylim, 'Color','c','LineWidth',2);

figure(2),hist(contrastB,50)
contrastBoot = mean(contrastB)
contrastSEBoot = sqrt(var(contrastB)) % no standard error except via Bootstrap
contrastBootBias = contrastBoot - contrast % huge!
contrastBiasCorrected = contrast - contrastBootBias % large correction
line([contrast contrast],ylim, 'Color','r','LineWidth',2);
line([contrastBoot contrastBoot],ylim, 'Color','m','LineWidth',2);
line([contrastBoot contrastBoot]-contrastSEBoot,ylim, 'Color','c','LineWidth',2);
line([contrastBoot contrastBoot]+contrastSEBoot,ylim, 'Color','c','LineWidth',2);




