% Exercise on comparison of GOF tests, NP problem 3.8
% fcp 140826, updated with in input 141030
%
%%
clear
set(0,'DefaultLineLineWidth',2);
set(0,'DefaultAxesFontSize',14);
set(0,'DefaultTextFontSize',18);

%%
mu = 0;
sigma = 1;
HWHM = sqrt(2*log(2));      % HWHM of standard normal
Nexpt = 100000;             % Number of experiments to simulate
Nbig = 100;                 % Sample size in an experiment
Nsmall = 10;                % sample size for small experiment

standardNormal = makedist('normal',0,1);

ad = zeros(1,Nexpt);
pH0ads = zeros(1,Nexpt);
ads = zeros(1,Nexpt);

pH0ks = zeros(1,Nexpt);
pH0ad = zeros(1,Nexpt);
pH1ks = zeros(1,Nexpt);
pH1ad = zeros(1,Nexpt);
pH1kss = zeros(1,Nexpt);
pH1ads = zeros(1,Nexpt);

adH1 = zeros(1,Nexpt);

% Accurate values of AD statistic
%adacc = zeros(1,Nexpt);
%adsacc = zeros(1,Nexpt);
%adH1acc = zeros(1,Nexpt);

%parfor nexpt = 1:Nexpt
for nexpt = 1:Nexpt
    unorm = normrnd(mu,sigma,Nbig,1);  
    ucauchy = HWHM*trnd(1,Nbig,1);     % Use t-distribution with 1 dof for Cauchy
%   first, try out tests on the null distribution sample
    [~,pH0ks(nexpt)] = kstest(unorm);
    % either adtest or adtestNP can be used here
    [~,pH0ad(nexpt),ad(nexpt)] = adtestNP(unorm,'Distribution',standardNormal);
    [~,pH0ads(nexpt),ads(nexpt)] = adtestNP(unorm(1:Nsmall),'Distribution',standardNormal);
    
%    adacc(nexpt) = adstatnormcdf(unorm);
%    adsacc(nexpt) = adstatnormcdf(unorm(1:Nsmall));
    
%   next, try the tests on the alternative sampling
    [~,pH1ks(nexpt)] = kstest(ucauchy);
    % must use adtestNP for sampling from this alternative to deal with
    % long tails
    [~,pH1ad(nexpt),adH1(nexpt)] = adtestNP(ucauchy,'Distribution',standardNormal);
    [~,pH1kss(nexpt)] = kstest(ucauchy(1:Nsmall));
    [~,pH1ads(nexpt)] = adtestNP(ucauchy(1:Nsmall),'Distribution',standardNormal);
    
%    adH1acc(nexpt) = adstatnormcdf(ucauchy);
end

%%
% Check behavior on H0 samples

nbins = 100;
figure;
h1 = subplot(1,2,1);
[histpH0ks] = hist(pH0ks,nbins);
hist(pH0ks,nbins);
xlabel('p value');
ylabel('Entries/0.01')
h = findobj(gca,'Type','patch');
set(h,'Facecolor',[0 .5 .5],'Edgecolor',[0 .5 .5]);

h2 = subplot(1,2,2);
[histpH0ad] = hist(pH0ad,nbins);
hist(pH0ad,nbins);
xlabel('p value');
h = findobj(gca,'Type','patch');
set(h,'Facecolor',[0 .5 .5],'Edgecolor',[0 .5 .5]);
axis([h1,h2],'square');

%%
% Do chi-square test for uniformity

expected = (Nexpt/nbins)*ones(1,nbins);

pchisqks = chi2cdf(sum((histpH0ks - expected).^2./expected),nbins-1,'upper')
pchisqad = chi2cdf(sum((histpH0ad - expected).^2./expected),nbins-1,'upper');


[pH0Sort,indexSort] = sort(pH0ad);
[pH0Sorts,indexSorts] = sort(pH0ads);     % for small smaple size
% for i = 1:Nexpt
%     adSort(i) = ad(indexSort(i));
%     adSorts(i) = ads(indexSorts(i));
% end
adSort = ad(indexSort);
adSorts = ads(indexSorts);

%adSortacc = adacc(indexSort);
%adSortsacc = adsacc(indexSorts);

%%
% critical value vs significance
figure;
semilogx(pH0Sort,adSort,'r',pH0Sorts,adSorts,'g');
xlim([0.0001,1]);
xlabel('Significance');
ylabel('Critical value');

%figure;
%semilogx(pH0Sort,adSortacc,'r',pH0Sorts,adSortsacc,'g');
%xlim([0.0001,1]);
%xlabel('Significance');
%ylabel('Critical value');
%title('Accurate estimates')

%%
% AD statistic under H1
figure;
xlo = 0;
dx = 5;
xhi = 2000;
bins = xlo:dx:xhi;
hist(adH1(adH1<xhi),bins);
xlim([0 xhi]);
xlabel('AD statistic');
ylabel('Entries/5');
h = findobj(gca,'Type','patch');
set(h,'Facecolor',[0 .5 .5],'Edgecolor',[0 .5 .5]);
minadH1 = min(adH1)
maxadH1 = max(adH1)

%figure
%hist(adH1acc(adH1acc<xhi),bins);
%xlim([0 xhi]);
%title('Accurate values')
%xlabel('AD statistic');
%ylabel('Entries/5');

%%
% Want to make a graph of power vs significance
% We only need the results from the H1-is-true studies,
% since we have the p-values (for H0).
% Let P be a given p value such that H0 will be rejected if p<P
% Count how many experiments have p<P; call this R(P) 
% Then power(P) = 1 - beta(P) = R(P)/Nexpt
% alpha = P is the significance

[ksecdf,kspv] = ecdf(pH1ks);
[adecdf,adpv] = ecdf(pH1ad);
[ksecdfs,kspvs] = ecdf(pH1kss);
[adecdfs,adpvs] = ecdf(pH1ads);

figure;
h2 = subplot(1,2,2);
plot(kspv,ksecdf,'b',adpv,adecdf,'r');
xlim([0,.1]);
xlabel('significance');

h1 = subplot(1,2,1);
plot(kspvs,ksecdfs,'b',adpvs,adecdfs,'r');
xlim([0,.1]);
xlabel('significance');
ylabel('power');

axis([h1,h2],'square');

figure;
h2 = subplot(1,2,2);
semilogx(kspv,ksecdf,'b',adpv,adecdf,'r');
xlim([0.0001,1]);
xlabel('significance');

h1 = subplot(1,2,1);
semilogx(kspvs,ksecdfs,'b',adpvs,adecdfs,'r');
xlim([0.0001,1]);
xlabel('significance');
ylabel('power');

axis([h1,h2],'square');

% Graph p-value against AD statistic 
idx = adH1(:) < 20;
figure;
semilogy(sort(adH1(idx)),sort(pH1ad(idx),'descend'));
xlabel('AD statistic');
ylabel('p value');

% Determine power at significance = 0.01
ksindex01 = find(kspv >= 0.01);
kssig = kspv(ksindex01(1));
kspower = ksecdf(ksindex01(1));

adindex01 = find(adpv >= 0.01);
if(size(adindex01) == [0 1])
    adsig = 0.01;
    adpower = 1;
else
    adsig = adpv(adindex01(1));
    adpower = ksecdf(adindex01(1));
end

ksindex01 = find(kspvs >= 0.01);
kssigs = kspvs(ksindex01(1));
kspowers = ksecdf(ksindex01(1));

adindex01 = find(adpvs >= 0.01);
adsigs = adpvs(adindex01(1));
adpowers = adecdfs(adindex01(1));

table(kssig, kspower, adsig, adpower, kssigs, kspowers, adsigs, adpowers)


