% Demonstrate calculation of confidence thresholds for pearson and kendall.
% THIS VERSION TAKES X DATA FROM A PUTATIVE 8-CHANNEL DEVICE, I.E. 8 DISCRETE DATA
% VALUES AVAILABLE FOR X, BUT FLOATING POINT DATA FOR Y. SO THERE ARE MANY TIES.
% For simplicity, look at just a single P-value. Note that the confidence
% is two-sided. 
function pval_ties1(Pvalue, Ntrials)

Confidence = 1 - Pvalue
Sigmas = erfinv(Confidence)*sqrt(2) % Number of sigmas equivalent for the normal distribution

% Begin with real, floating point data drawn from a uniform distribution. 
% We expect no "ties", even with large datasets.
% The distribution does not matter at all for uncorrelated data (right?).

Ndata = 8:2:30;
Npopulationsizes = numel(Ndata); % number of different population sizes to test via MonteCarlo
Ntrials % number of trials to run for each population size
for k=1:Npopulationsizes
   PopulationN = Ndata(k) % diagnostic
   tic % diagnostic
   for trial = 1:Ntrials
      x = floor(8*rand(Ndata(k),1));
      y = rand(Ndata(k),1);
      rho(trial) = spearman(x,y);
      tau(trial) = kendall(x,y);
   end
   rho = sort(abs(rho));
   tau = sort(abs(tau));
   rho_threshold(k) = rho( ceil(Confidence*Ntrials) );
   tau_threshold(k) = tau( ceil(Confidence*Ntrials) );
   rho_threshold_normal(k) = Sigmas*sqrt( 1 / (Ndata(k)-1) ); 
      % per variance given by Octave: help spearman
   tau_threshold_normal(k) = Sigmas*sqrt( (2 * (2*Ndata(k)+5)) / (9 * Ndata(k) * (Ndata(k)-1)) ); 
      % per variance given by Octave: help kendall
   toc % diagnostic
end


Confidence_string = sprintf("%5.2f", Confidence*100)
Ntrials_string = sprintf("%d", Ntrials)

figure(1)
hold off
plot( Ndata, rho_threshold, 'b-o', Ndata, rho_threshold_normal, 'k-+')
hold on
title(['8-channel data: Spearman \rho, ', Confidence_string, '% Significance, ', Ntrials_string, ' Trials'])
xlabel('Population Size (N)')
ylabel('|\rho| threshold')
legend('MonteCarlo','Normal Assumption')
print(['Spearman-ties1-', Confidence_string, 'conf-', Ntrials_string, 'trials.pdf'], '-dpdfwrite')

figure(2)
hold off
plot( Ndata, tau_threshold, 'b-o', Ndata, tau_threshold_normal, 'k-+')
hold on
title(['8-channel data: Kendall \tau, ', Confidence_string, '% Significance, ', Ntrials_string, ' Trials'])
xlabel('Population Size (N)')
ylabel('|\tau| threshold')
legend('MonteCarlo','Normal Assumption')
print(['Kendall-ties1-', Confidence_string, 'conf-', Ntrials_string, 'trials.pdf'], '-dpdfwrite')
