%% Hypothesis Testing
%
% Python notebook for illustrating the concept of Hypothesis Testing and specific test statistics; among them the very useful Kolmogorov-Smirnov test.
% The Kolmogorov-Smirnov test (KS-test) is a general test to evaluate if two distributions in 1D are the same. This program applies an unbinned KS test, and compares it to a $\chi^2$-test and a simple comparison of means. The distributions compared are two unit Gaussians, where one is then modified by changing:
% Mean
% Width
% Normalisation
% The sensitivity of each test is then considered for each of these changes.
%% References:
% Barlow: p. 155-156
% __[Wikipedia: Kolmogorov-Smirnov test](http://en.wikipedia.org/wiki/Kolmogorov-Smirnov_test)__
% Though influenced by biostatistics, a good discussion of p-values and their distribution can be found here:
% [How to interpret a p-value histogram?](http://varianceexplained.org/statistics/interpreting-pvalue-histogram/)
%
%% Authors:
% Troels C. Petersen (Niels Bohr Institute)
%% Date:
% 03-12-2020 (latest update)
%

N_exp = 1;
N_events_A = 100;
N_events_B = 100;

dist_mean_A  =  0.0;
dist_width_A =  1.0;
dist_mean_B  =  0.0;
dist_width_B =  1.0;

N_bins = 100;
xmin = -5.0;
xmax = 5.0;

all_p_mean = zeros([1,N_exp]);
all_p_chi2 = zeros([1,N_exp]);
all_p_ks   = zeros([1,N_exp]);

verbose=true;

for iexp=1:N_exp
%     if ((iexp+1)%1000 == 0):
%         print(f"Got to experiment number: {iexp+1}")
%     end
    % Generate data:
    x_A_array = normrnd(dist_mean_A, dist_width_A, [1, N_events_A]);
    x_B_array = normrnd(dist_mean_B, dist_width_B, [1, N_events_B]);
    
    
    % Test if there is a difference in the mean:
    % ------------------------------------------
    % Calculate mean and error on mean:
    [mean_A, width_A, sdom_A] = mean_std_sdom(x_A_array);
    [mean_B, width_B, sdom_B] = mean_std_sdom(x_B_array); 

    % Consider the difference between means in terms of the uncertainty:
    d_mean = mean_A - mean_B;
    % ... how many sigmas is that away?

    % Turn a number of sigmas into a probability (i.e. p-value):
    p_mean    = 0.5;      % Calculate yourself. HINT: "stats.norm.cdf or stats.norm.sf may be useful!"
    all_p_mean(iexp) = p_mean;
    
    
    % Test if there is a difference with the chi2:
    % --------------------------------------------
    % Chi2 Test:
    p_chi2 = 0.5;         % Calculate the p-value of the Chi2 between histograms of A and B yourself.
    all_p_chi2(iexp) = p_chi2;

    
    % Test if there is a difference with the Kolmogorov-Smirnov test on arrays (i.e. unbinned):
    % -----------------------------------------------------------------------------------------
    [h, p_ks] = kstest2(x_A_array, x_B_array);           % Fortunately, the K-S test is implemented in stats!
    all_p_ks(iexp) = p_ks;


    % Print the results for the first 10 experiments
    if (verbose) && (iexp < 10)
      disp(num2str(iexp)+":  p_mean: " +num2str(p_mean)+  ", p_chi2: "+num2str(p_chi2)+", p_ks: "+num2str(p_ks))
    end
    % In case one wants to plot the distribution for visual inspection:
    if (iexp == 1)
        edges = [xmin:(xmax-xmin)/N_bins:xmax];
        histogram(x_A_array, edges)%, histtype='step', label='A', color='blue')
        %ax1.set(title='Histograms of A and B', xlabel='A / B', ylabel='Frequency / 0.05')  
        hold on
        histogram(x_B_array, edges)%, histtype='step', label='B', color='red')

    end
end

%

N_bins = 50

if (N_exp > 1)
    figure()
    fontSize = 10;
    edges = [0:1/N_bins:1];
    histogram(all_p_mean, edges)
    title('Histogram, probability mu', 'FontSize', fontSize);
    xlabel('p-value', 'FontSize', fontSize);
    ylabel('Frequency / 0.02', 'FontSize', fontSize);
    figure()
    histogram(all_p_chi2, edges)
    title('Histogram, probability chi2', 'FontSize', fontSize);
    xlabel('p-value', 'FontSize', fontSize);
    ylabel('Frequency / 0.02', 'FontSize', fontSize);
    figure()
    histogram(all_p_ks, edges)
    title('Histogram, probability Kolmogorov', 'FontSize', fontSize);
    xlabel('p-value', 'FontSize', fontSize);
    ylabel('Frequency / 0.02', 'FontSize', fontSize);
end

%
%
% Define a function that calculates the mean, standard deviation and the standard deviation (i.e. uncertainty) on mean (sdom):

function [mean_x, std_x, std_sqrt]=mean_std_sdom(x)
    std_x = std(x);
    mean_x = mean(x);
    std_sqrt = std_x/sqrt(length(x));
end 

%%  Questions:
%
% 1. First run the program with one experiment (N_exp = 1) to display the two distributions A and B, when:
%     - They are the same.
%     - The mean of A is increased (to e.g. 0.1).
%     - The width of A is enlarged (to e.g. 1.2).
%     - The normalisation of A is increased.
%
% Get a feel for how much you need to change the distribution, before you can _by eye_ see that they are not the same. I.e. can you see any difference, if `mean_A` $= 0.1$? Or how about $0.2$? How do you quantify this and when do you start to doubt? And how about `width_A` $= 1.1$? Or $1.2$? Again, can you see it by eye? Finally, try to put $1050$ events into B. Is that visible? How about $1100$? And do you see an impact in the p-values?
%
% 2. Could you for the test of the means have calculated how much of a change in the mean is needed for a difference to be statistically significant? Do so, and see if it somewhat matches you guess/estimate from above!
%
%
% 3. Now run the tests 1000 times, where A and B are unit Gaussians and thus identical. How should the distributions of the test probabilities come out? And is this the case, approximately? If not, think of reasons for this, and what could be a remedy. HINT: Large statistics is always easier!
%
%
% 4. Repeat the changes in question 1), and see which tests "reacts" most to these modifications. How much of a change in the mean is required for 95% of the tests (of each kind) to give a probability below 5%? How much is required for the width? And the norm?
%
%
% 5. Possibly try to test different distributions than the Gaussian one (e.g. exponential, uniform, etc.), and see how the tests performs.
%
%
% NOTE: The Kolmogorov-Smirnov test has the great advantage that it can handle ANY distribution (even the Cauchy distribution - remind yourself of that one!). The reason is, that it doesn't care about any PDF, nor how far out an outlier is. It is just a matter of the difference in integrals between the two functions.
%
%
%% Advanced:
%
% 6. Obviously, the test of the means is not sensitive the a change in the width. Make such a test yourself by calculating the widths and the uncertainty on the widths (or perhaps try the F-test!). Note that in a (unit) Gaussian the uncertainty on the width is of the same order as that of the means!
%
%
%% Very advanced:
% 7. Implement in python the following tests:
%      - Lilliefors test
%      - Shapiro-Wilk test
%      - Anderson-Darling test
%      - Cramer-von-Mises test
%      - Jarque-Bera test
%      - Kuiper's test
%      - Mann-Whitney-Wilcoxon test
%      - Siegel-Tukey test
%
% and quantify under various conditions and datasets the power of each and the correlation among them. Write it up, and send it to a statistics journal. :-)

