% quick_ttest_RP26Apr2016.m
%
% Simple program to run t-tests on simulated non-Gaussian data, for
% illustrating problems with p-values, NHST etc.
% Simulated datasets: log-normal distributed numbers, i.e. 
%    exp(Gaussian distributed numbers), i.e. exp(N(m, s^2)).
%
% Inputs
%   N : number of points in each trial (e.g. 10)
%   Ntrials : number of trials
%   m0, m1 : mean of Gaussian distribution for set 0, 1
%   std0, std1 : standard deviation of Gaussian distribution for set 0, 1
%   plotopt : if true, make a graph
%
% Outputs
%   e0, e1 : simulated data points, N x Ntrials matrix
%   greatermean : Ntrial x 1 boolean array, true if mean(e1) > mean(e0)
%   p = p-value from a t-test, Ntrial x 1 array
%   p_ks = p-value from a Kolmogorov-Smirnoff test, Ntrial x 1 array
%
% Raghuveer Parthasarathy
% Apr 26, 2016
% last modified Sept. 22, 2016


function [e0, e1, greatermean, p, p_ks] = ...
    quick_ttest_RP26Apr2016(N, Ntrials, m0, std0, m1, std1, plotopt)

rng('shuffle'); % initialize random numbers

%% Defaults
if ~exist('N', 'var') || isempty(N)
    N = 10;
end
if ~exist('Ntrials', 'var') || isempty(Ntrials)
    Ntrials = 100;
end
if ~exist('m0', 'var') || isempty(m0)
    m0 = -0.1; % mean of distr. 0, before exp.
end
if ~exist('m1', 'var') || isempty(m1)
    m1 = 0.4; % mean of distr. 1, before exp.
end
if ~exist('std0', 'var') || isempty(std0)
    std0 = 0.5; % std. dev. of distr. 0, before exp.
end
if ~exist('std1', 'var') || isempty(std1)
    std1 = 1.5; % std. dev. of distr. 1, before exp.
end
if ~exist('plotopt', 'var') || isempty(plotopt)
    plotopt = true;
end

%% Draw all the random numbers
e0 = exp(m0 + std0*randn(N,Ntrials));
e1 = exp(m1 + std1*randn(N,Ntrials));

% Calculate all the p-values, etc.
greatermean = zeros(Ntrials,1);
p = zeros(Ntrials,1);
p_ks = zeros(Ntrials,1);
for j=1:Ntrials
    greatermean(j) = mean(e1(:,j)) > mean(e0(:,j)); % true if mean(e1) > mean(e0)
    [~,p(j)] = ttest2(e0(:,j), e1(:,j)); % t-test
    [~,p_ks(j)] = kstest2(e0(:,j), e1(:,j)); %Kolmogorov-Smirnov test
end

% Just look at the first draw
disp('One trial...')
fs = sprintf('  p = : %.03f', p(1)); disp(fs);
fs = sprintf('  Is mean_1 > mean_0? : %d', greatermean(1)); disp(fs);

if plotopt
    figure; plot(0.2*(rand(N,1)-0.5) + ones(N,1), e0(:,1), 'ko', ...
        'markerfacecolor', 0.7*[1 1 1], 'markersize', 10)
    hold on
    plot(0.2*(rand(N,1)-0.5) + 2*ones(N,1), e1(:,1), 'bd', ...
        'markerfacecolor', [0.2 0.5 1], 'markersize', 10)
    a = axis;
    axis([0 3 0 a(4)])
    ylabel('Measurement')
    set(gca, 'XTickLabel', '')
    %text(0.75, 1, 'Treatment 1', 'rotation', 90)
    %text(1.75, 1, 'Treatment 2', 'rotation', 90)
    text(0.6, -0.7, 'Treatment 1')
    text(1.6, -0.7, 'Treatment 2')
end

outputcsv = false;
if outputcsv
    e_both = zeros(N, 2*Ntrials);
    e_both(:,1:2:end) = e0;
    e_both(:,2:2:end) = e1;
    csvwrite('ttest_RP_measurements.csv', e_both);
end

    
if Ntrials>1
    figure('position', [600 100 500 500]); histogram(p)
    xlabel('p value')
    ylabel('Counts')
    fs = sprintf('Number of trials with mean_1 > mean_0: %d out of %d', ...
        sum(greatermean==1), Ntrials);disp(fs)
    fs = sprintf('Number of trials with p < 0.05: %d out of %d', ...
        sum(p<0.05), Ntrials); disp(fs);
    fs = sprintf('Number of trials with p_ks < 0.05: %d out of %d', ...
        sum(p_ks<0.05), Ntrials); disp(fs)
end

end
