%% ChiSquare Test:
% This Python macro illustrates the use of ChiSquare as a goodness-of-fit measure,
% how the ChiSquare distribution comes about, how to use the ChiSquare value (along with
% number of degrees of freedom, Ndof) as a test, and that it actually works!
% The program generates a certain number of datasets, each consisting of 9 points along
% a line. These are then fitted, and the result and the Chi2 of the fit is recorded along
% with the probability of the fit.
%% References:
% - Barlow: Chapter 6
% - Cowan: Chapter 2.7, Chapter 7
% - Bevington: Chapter 6
%% Author(s), contact(s), and dates:
% - Author: Troels C. Petersen (NBI)
% - Email:  petersen@nbi.dk
% - Date:   12th of November 2020
%
%% Initial Question/Task:
% Make sure you've read the relevant references and that you understand not only what the ChiSquare is, but also that it follows the ChiSquare distribution, and that the probability of obtaining such a ChiSquare or worse (assuming that the hypothesis/PDF is the correct one!) can be calculated from it.
% Executive Summary: The ChiSquare method is used both for getting a good and robust fit, and for testing the quality/probability of this fit.
%


Nexp = 1000;
Npoints = 9;
alpha0 = 3.6;
alpha1 = 0.3;
alpha2 = 0.0;    % Initially, we don't include a quadratic term!
sigmay = 1.0;     % The uncertainty on each point is 1.0


% Arrays for storing fit results:
array_alpha0 = zeros(Nexp);
array_alpha1 = zeros(Nexp);
array_Chi2 = zeros(Nexp);
array_Prob = zeros(Nexp);

%# Define a fit function:

%


array_alpha0 = [];
array_alpha1 = [];
array_Chi2 = [];
array_Prob = [];
%Loop, repeating the data generation and fit:
for iexp=1:Nexp
    
    % Generating data by filling values into (x,y) and associated uncertainties (here chosen to be 0 for x):
    x = 1:Npoints+1;
    ex = zeros(x);
    y = alpha0 + alpha1*x + normrnd(0,sigmay,[1,Npoints+1]) + alpha2*x.^2;
    % Note how we include uncertainty in y - by simply adding a Guassian number to the "theoretical value".
    ey = sigmay*ones(1,length(x));

    %# ------------------------------------------------------------------ #
    %# Numerical fit:
    %# ------------------------------------------------------------------ #

    func = @(a) a(1) + a(2)*x;
    chi2 = @(a) sum(((y-(func(a)))./ey).^2);
    %# Here we let Minuit know, what to minimise, how, and with what starting parameters:   
    guess = [alpha0, alpha1];
    [optim,chi_2,exitflag,output] = fminsearch(chi2,guess);
      
    Nvar = 2;                     %# Number of variables (alpha0 and alpha1)
    Ndof_fit = Npoints - Nvar;    %# Number of degrees of freedom = Number of data points - Number of variables
    prob = 1-chi2cdf(chi_2,Ndof_fit); % probability value 1-cdf
    
    array_alpha0(iexp) = optim(1);
    array_alpha1(iexp) = optim(2);
    array_Chi2(iexp) = chi_2;
    array_Prob(iexp) = prob;
end

%
%% Fit the data and plot the result:
% Below we plot the latest fit (i.e. the last one from above), so that we can actually see what is going on (and potentially wrong!).

errorbar(x,y, ey, '.')
hold on; plot(x, func(optim))
xlabel('x'); ylabel('y');
disp("alpha1: " + num2str(optim(1)) +newline+ ...
    'alpha2: ' + num2str(optim(2)) + newline+ ...
    'Chi2: ' + num2str(chi_2) + newline+ ...
    'Ndof: ' + num2str(Ndof_fit) + newline+ ...
    'Prob: ' + num2str(prob));

%
%% Many experiments
% In case we have more than one "experiment" (i.e. fit), we would like to see how the fit results and the Chi2 values distribute themselves:

%# Distribution of alpha0 values:
edges = [0:6/60:6];
figure; histogram(array_alpha0, edges)
%# Distribution of alpha1 values:
edges = [-0.2:0.8/50:0.8];
figure; histogram(array_alpha1)%
%# Distribution (normed) of Chi2 values, along with Chi2-distribution (for Ndof=7) drawn on top:
edges = [0:25/25:25.0];
figure; histogram(array_Chi2)
%# Distribution of Chi2 fit probabilities:
edges = [0:1/20:1];
figure; histogram(array_Prob)

%% Remark on $\chi^2$ rule of thumb:
% A general rule is that the "reduced chisquare" (defined as $\chi^2$/Ndof) should be "around 1". That is true, but the rule definitely has its limitations, since "around" depends on the size of Ndof! This is examplified below, considering two cases with the same reduced chisquare of 1.5, but with **very** different conclusions.
% The reason for this can be found in the mean and width of the $\chi^2$ distribution (with Ndof = k), which are k and sqrt(2k). From the central limit theorem, the $\chi^2$ should approach a Gaussian for "large" values of k (k > 50 is probably enough), but the **relative** width (i.e. width/mean) becomes smaller and smaller as k grows, and a reduced chisquare of e.g. 1.5 becomes more and more unlikely, as k increases.

%# Numbers to illustrate the point, that using the ratio Chi2/Ndof ("reduced chisquare") is only a rule of thumb:
disp("Considering the relation between Chi2/Ndof ratio and Prob(Chi2,Ndof) vs. Ndof:"+newline+ ...
    "Low number of degrees of freedom:   Prob(  3.0,  2) = "+num2str(1-chi2cdf(3.0,2))+newline+ ...
    "High number of degrees of freedom:  Prob(300.0,200) = "+num2str(1-chi2cdf(300.0, 200)))

% The conclusion is, that while the reduced $\chi^2$ is a good rule of thumb, it is NOT the correct way to consider a $\chi^2$ value.
%
% ----------------------------------------------------------------------------------
%% Questions:
%  1) Run the code such that it does exactly one fit (initial setting), and take a look at the line fitted.
%     Does this look reasonable, and are you satisfied with the fit probability? Does the fit reproduce the input
%     numbers well (include the uncertainties in your answer)?
%
%  2) What is the chance that the input for the data could actually be from a flat distribution, i.e. without the slope?
%     [There is a correct (slower) way of doing this, and another quick but approximate way of testing this].
%
%  3) Now increase the number of experiments to e.g. 1000, and rerun the macro. Figure
%     out what you see in the window split in 2-by-2, and go through each of these to
%     see, if you understand every feature of the distributions shown, and if you are
%     happy with them! Is this what you expected? Are they also as wide as you would expect?
%     Possibly try to make some fits of these output distributions to check, if you're satisfied!
%     <br>
%     This somehow makes this the "long" question without any leading questions.
%     But you should try to be statistically minded enough and have an idea of what to
%     look for, at least to some degree. We will discuss this further in class.
%
%  4) Investigate if the distributions of probabilities is flat, or if it has some
%     slope to it. Do you understand why the distributions of probabilities is flat?
%     (This is typically conceptually hard to begin with, but don't worry!)
%
%  5) Find the line of code where the random points for the line are generated, and add a
%     quadratic term in x with the coefficient alpha2 = -0.04. Run the program again with Nexp=1,
%     and start by looking at the single fit in the graph. Can you see this change?
%     Now run 1000 experiments again, and see what has changed in the distributions
%     in the 2-by-2 window when including such a term. Are these changes as you would
%     have expected them?
%     <br>
%     NOTE: The reason why the quadratic term is more clear in the latter case is NOT
%     that it is a better method. It is simply because there is 1000 times more statistics!
%
%  6) With the quadratic coefficient from question 5), the linear fit of course does
%     not do very well. What changes are needed for the fit to be good again? Make
%     these changes, and see that the condition in question 4) is again met.
%
%%     (Semi-)Advanced questions:
%
%  7) Try to generate data according to more advanced functions, and then fit the data with these.
%     Do you still obtain a flat distribution in probability? Do you also manage to fit these with
%     simpler/other functions and still maintain the flat distribution?
%

