%power1.m
%Monte Carlo power analysis as a function of effect size
%and joystick time quantization sampling effects
%Power appears to be extremely robust across time quantization

clear all;

%define sampling time quantization error
time_quantum = 1;

%pick a group sample size
samplesize = 22;

%define theoretical treatment and control group populations
mean_treatment = 880;
sd_treatment = 141;
mean_control = 1018;  %very large effect size
%mean_control = 920;  %small effect size
sd_control = 192;

%seed random number generator
if verLessThan('matlab', '8.1.0')
    %RandStream.setDefaultStream(RandStream('mt19937ar','seed',sum(clock*100)));
    RandStream.setDefaultStream(RandStream('mt19937ar','seed',0));
else
    %RandStream.setGlobalStream(RandStream('mt19937ar','seed',sum(clock*100)));
    RandStream.setGlobalStream(RandStream('mt19937ar','seed',0));
end

%generate bootstrap distribution for t when null hypothesis is true
n_iterations = 20000;
bs_null_distribution = zeros(n_iterations,1);
for count = 1:n_iterations
   %draw two independent samples from control group's model distribution
   %this simulates conditions when the null hypothesis is true
   treatment_sample = normrnd(mean_control, sd_control, samplesize, 1);
   control_sample = normrnd(mean_control, sd_control, samplesize, 1);
   %apply usb sampling time quantization error to sample data
   s1 = ceil(treatment_sample/time_quantum) * time_quantum;
   s2 = ceil(control_sample/time_quantum) * time_quantum;
   %compute difference between sample means; add into bootstrap distribution
   bs_null_distribution(count) = mean(s1)- mean(s2);
end

%generate bootstrap distribution for t when null hypothesis is false
bs_alt_distribution = zeros(n_iterations,1);
for count = 1:n_iterations
   %draw an independent sample from each group's model distribution
   %this simulates conditions when the null hypothesis is false
   treatment_sample = normrnd(mean_treatment, sd_treatment, samplesize, 1);
   control_sample = normrnd(mean_control, sd_control, samplesize, 1);
   %apply usb sampling time quantization error to sample data
   s1 = ceil(treatment_sample/time_quantum) * time_quantum;
   s2 = ceil(control_sample/time_quantum) * time_quantum;
   %compute difference between sample means; add into bootstrap distribution
   bs_alt_distribution(count) = mean(s1)- mean(s2);

end

%find cutoff value for lower 5% of null hypothesis bootstrap distribution
[nullcumprob, nullxval] = ecdf(bs_null_distribution);
temp = abs(nullcumprob-0.05);
index=find(temp==min(temp));
nullcritx = nullxval(index(1));

%find beta probability for alt. hyp. bootstrap distribution
[altcumprob, altxval] = ecdf(bs_alt_distribution);
temp = abs(altxval-nullcritx);
index = find(temp==min(temp));
beta = 1-altcumprob(index(1));

%plot bootstrap distributions
%first, plot null hypothesis distribution
figure(1), clf
subplot(2,1,1);
hist(bs_null_distribution,33); %33 intervals
hold on;
h = findobj(gca, 'Type', 'patch');
set(h, 'FaceColor', 'w', 'EdgeColor', 'b');
axis([-350 200 0 4000]);
xlabel('Treatment-Control Mean (msec)');
ylabel('Frequency');
title(['Null Hypothesis, N = ', num2str(samplesize),', Quantization = ',num2str(time_quantum)]);
%annotate figure
f21str(1)={['Critical Mt-Mc (p < 0.05) = ',num2str(nullcritx)]};
text(-325,3500,f21str);
line([nullcritx nullcritx],[0 4000]);
hold off;
%
%next, plot the alternative hypothesis distribution
subplot(2,1,2);
hist(bs_alt_distribution,33); %33 intervals
hold on;
h = findobj(gca, 'Type', 'patch');
set(h, 'FaceColor', 'w', 'EdgeColor', 'b');
axis([-350 200 0 4000]);
xlabel('Treatment-Control Mean (msec)');
ylabel('Frequency');
title(['Alt. Hypothesis, N = ', num2str(samplesize),', Quantization = ',num2str(time_quantum)]);
effectSize=(mean_treatment-mean_control)/((sd_treatment+sd_control)/2);
%annotate figure
f22str(1)={['Mt-Mc = ',num2str(mean_treatment-mean_control),' msec']};
f22str(2)={['Effect Size = ',num2str(effectSize)]};
f22str(3)={['Beta = ',num2str(beta)]};
f22str(4)={['Power = ',num2str(1-beta)]};
f22str(5)={['Total N = ',num2str(samplesize*2)]};
text(0,2500,f22str);
line([nullcritx nullcritx],[0 4000]);
hold off;
