%% Hypothesis test of normally distrbuted data.
% A small study of various methods for testing whether observed data are compatible with a normal distribution, without knowledge of the distribution mean and standard deviation apart from what can be learned from the data temselves.
% Start by testing a dataset colleced by a student in 2019. The measurements are of the period of a pendulum with a stopwatch.

% Some student stopwatch data from 2019, measurements og period of pendulum
% in seconds.
t=[1.30 1.36*ones(1,3) 1.38*ones(1,6) 1.40*ones(1,6) 1.42*ones(1,6) 1.44*ones(1,3) 1.46 1.48 1.48 1.50 1.52];

% Normalize the data to the mean of the data and the standard deviation of
% the data.
tnorm=(t-mean(t))/std(t);
histogram(x,10)
fprintf('Mean period = %.3f s\n',mean(t))
fprintf('Std. deviation of the periods = %.3f s\n',std(t))
fprintf('Uncertainty on the mean period = %.3f s\n',std(t)/sqrt(length(t)));

%
%% kstest on single observed distribution
% Perform the Kolmogorov-Smirov test, which tests whether observed data are compatible with a unit normal distribution (the null hypothesis). Since these data are unnormalized the mean is not 0 and the standard deviation is not 1, so the test should result in rejecting the null hypothesis at the specified confidence level (95% by default).

% If the data are significantly (95% CL) NOT compatible with a unit normal
% distribution(95% CL) the kstest will return False (0)
if kstest(t)
    fprintf('The UNnormalized data appear to NOT be unit-normally distributed at the 95%% CL.\n')
else
    fprintf('The UNnormalized data appear to be unit-normally distributed at the 95%% CL.\n')
end

%
% Repeat the test on the data normalized to the observed mean and observed standard deviation.

% If the data are significantly (95% CL) NOT compatible with a unit normal
% distribution (95% CL) the kstest will return False (0)
if kstest(tnorm)
    fprintf('The normalized data appear to NOT be unit-normally distributed at the 95%% CL.\n')
else
    fprintf('The normalized data appear to be unit-normally distributed at the 95%% CL.\n')
end

%
%% Sample properties of kstest
% So far we have only tested a single dataset. If the confidence level is to have any meaning then in an ensemble of experiments (or simulated experiments, also referred to as pseudo-experiments) we expect the null hypotheis to be rejected incorretly at a rate of 1-CL, so 5% by default. The following simulation shows that this is not the case when the data-driven estimates of the mean and standard deviation are used to normalize the data, but the 1-CL rejection is acheived when the data are normalized to the known mean and standard deviation. There turns out to be no correspondence between the confidence level and the rejection rate. This also means that this method of detecting defviations from unit normality is not sensitive.

% Simulate a distribution of measured periods with the mean and standard
% deviation from a student's 2019 data.

m = 1.411; % mean of periods
s = 0.046; % standard deviation of periods
N = 30;    % Number of period measurements
M = 10000;    % Number of pseudo-experiments (simulated experiments)

% Perform M pseudo-experiment and see how often kstest of 
%  * the normalized data
%  * the unnormalized data
%  * the data normalized to the known mean and standard deviation
% fails.

normRej = 0;
unnormRej = 0;
normRej0 = 0;
pnorm = zeros(M,1);
pnorm0 = zeros(M,1);
CL = 0.95;
for i=1:M
    t = normrnd(m,s,N,1);
    tnorm0=(t-m)/s;
    tnorm=(t-mean(t))/std(t);
    [hnorm,pnorm(i)] = kstest(tnorm);
    [hnorm0,pnorm0(i)] = kstest(tnorm0);
    if kstest(t)    
        unnormRej = unnormRej + 1;
    end
    if hnorm0
        normRej0 = normRej0 + 1;
    end
    if hnorm
        normRej = normRej + 1;
    end
end

%
% The two histograms show (1) that the distribution of p-values using the data normalized to themselves is skewed toward high values and (2) that the distribution of p-values normalized to the true mean and standard deviaiton is approximately uniform. In other words, the kstest on the data normal.ized to themselves will reject the unit-normality hypothesis much less frequently than the intended 1-CL rate given a true hypothesis and will fail to reject the unit-normality hypothesis of non-unit-normal distributions at a high rate.

histogram(pnorm)
title('Pseudo-data normalized to themselves')
xlabel('p-value'), ylabel('Entries')
histogram(pnorm0)
title('Pseudo-data normalized to true values')
xlabel('p-value'), ylabel('Entries')
fprintf('Rate of rejecting UNnormalized data: %.1f %%\n',100*unnormRej/M);
fprintf('Rate of rejecting normalized   data: %.1f %%\n',100*normRej/M)
f = normRej0/M;
df = sqrt(f*(1-f)/M);
fprintf('Rate of rejecting data normalized to true mean and std.dev.: %.1f +- %.1f %%\n',100*f,100*df)
fprintf('Note that 1-CL is %.1f %%\n',100*(1-CL))

%% Lilliefors hypothesis test
% Two other tests are more appropriate for testing data with unknown mean and standard deviation: The Andersen-Darling goodness-of-fit hypothesis test (adtest(x), where x are the raw data), and the Lilliefors goodnes-of-fit test of composite normality (lillietest(x), where x are the raw data).
% The following simulation shows that the Lilliefors test gives approximately the expected rejection rate of the null hypothesis, regardless of whether the data are normalized or not.
%

normRej = 0;
unnormRej = 0;
normRej0 = 0;
for i=1:M
    t = normrnd(m,s,N,1);
    tnorm0=(t-m)/s;
    tnorm=(t-mean(t))/std(t);
    if lillietest(t)
        unnormRej = unnormRej + 1;
    end
    if lillietest(tnorm0)
        normRej0 = normRej0 + 1;
    end
    if lillietest(tnorm)
        normRej = normRej + 1;
    end
end
f = unnormRej/M;
df = sqrt(f*(1-f)/M);
fprintf('Rate of rejecting UNnormalized data: %.1f +- %.1f %%\n',100*f,100*df);
fprintf('Rate of rejecting data: %.1f %%\n',100*normRej/M)
fprintf('Rate of rejecting data normalized to observed mean and std.dev.: %.1f %%\n',100*normRej0/M)

%
%% Conclusion
% Don't use the Kolmogorov-Smirnov test (kstest) unless you want to test a specific hypothesis with a known mean and standard deviation (and you must normalize the data to these parameters). If you nornalize the data to the observed mean and standard deviation the kstest will be less powerful than Andersen-Darling or Lilliefors tests. These latter two tests take as input raw data and give results that respect the stated confidence level of the hypothesis test.