% enee324StatsPermTestEx.m

close all; clear all
% can get file via web browser, or via this command
% urlwrite('http://terpconnect.umd.edu/~jzsimon/enee324/hw/permExData.mat','permExData.mat')
load('permExData.mat')

% This example and the data come from
% Bootstrap Methods and Permutation Tests by Hesterberg, Moore, Monaghan, Clipson, and Epstein
% http://home.comcast.net/~timhesterberg/bootstrap/
%
% TELECOMMUNICATION REPAIR TIMES
%
% Verizon is the primary local telephone company (the legal term is Incumbent Local
% Exchange Carrier, ILEC) for a large area in the eastern United States. As such,
% it is responsible for providing repair service % for the customers of other 
% telephone companies (known as Competing Local % Exchange Carriers, CLECs) in this
% region. Verizon is subject to fines if the repair times (the time it takes to fix
% a problem) for CLEC customers are substantially worse than those for Verizon?s own
% customers. This is determined using hypothesis tests, negotiated with the local 
% Public % Utilities Commission (PUC).

data1 = ILEC_data;
data2 = CLEC_data;
n1 = length(data1)
n2 = length(data2)

figure % explore data 1 probability distribution
hist(data1,20) % not enough bins, so
hist(data1,100)

figure % explore data 2 probability distribution
hist(data2,20) % too many bins, so
hist(data2,10)

% put both probability distribution figures on equal footing
figure(1);xlim([0 200])
figure(2);xlim([0 200])

% calculate means and difference
m1 = mean(data1)
m2 = mean(data2)
dm = m2-m1

% calculate p-value via t-test
[h,p_ttest]=ttest2(data1, data2, 0.01, 'left');
p_ttest

%Permutation data framework
datapool = [data1;data2];
nTotal = length(datapool);
nPT = 5000-1; % so nPT + 1 is a nice round number

data1PT = nan(size(data1)); 
data2PT = nan(size(data2));
dmPT = nan(1,nPT);

for k = 1:nPT
    datapoolPT = datapool(randperm(length(datapool)));
    data1PT = datapoolPT(1:n1);
    data2PT = datapoolPT(n1+1:end);
    
    m1PT = mean(data1PT);
    m2PT = mean(data2PT);
    dmPT(k) = m2PT-m1PT;
end
p_PT_dm = 1-(sum(dm>=dmPT)+1)/(nPT+1) % calculate PT p-value

% done, but see results visually too
dmsorted = sort(dmPT);
dm95=dmsorted(0.95*(nPT+1));
dm99=dmsorted(0.99*(nPT+1));

dmtt=dmsorted(round((1-p_ttest)*(nPT+1))); % same for t-test result

figure;hist(dmPT,50);title('mean difference permutation distribution')
 line(dm95*[1 1],[0 nPT/10],'linewidth',1,'color',[1 0 1])
 line(dm99*[1 1],[0 nPT/10],'linewidth',1,'color',[1 0 0])
 line(dm*[1 1],[0 nPT/10],'linewidth',1,'color',[0 1 0])
 line(dmtt*[1 1],[0 nPT/10],'linewidth',1,'color',[0 1 1])
 legend('permutation values','95% significance','99% significance','measured value',...
     't-test equivalent (incorrect)',0)
 
% optional: explore Permutation CDF
figure;
subplot(2,1,1)
plot(dmsorted,[1:nPT]/nPT);title('mean difference permutation CDF')
subplot(2,1,2)
plot(dmsorted,[1:nPT]/nPT);title('Expanded mean difference permutation CDF')
 ylim([0.90 1]);
 line(dm95*[0 1],0.95*[1 1],'linewidth',1,'color',[1 0 1])
 line(dm99*[0 1],0.99*[1 1],'linewidth',1,'color',[1 0 0])
 line(dm*[0 1],(1-p_PT_dm)*[1 1],'linewidth',1,'color',[0 1 0])
 line(dmtt*[0 1],(1-p_ttest)*[1 1],'linewidth',1,'color',[0 1 1])
 line(dm95*[1 1],0.95*[0 1],'linewidth',1,'color',[1 0 1])
 line(dm99*[1 1],0.99*[0 1],'linewidth',1,'color',[1 0 0])
 line(dm*[1 1],(1-p_PT_dm)*[0 1],'linewidth',1,'color',[0 1 0])
 line(dmtt*[1 1],(1-p_ttest)*[0 1],'linewidth',1,'color',[0 1 1])
 legend('permutation values','95% significance','99% significance','measured value',...
     't-test equivalent (incorrect)',0)


