% Paired Samples - Matlab t-test for Paired Samples
% Barukyah Shaparenko and John M. Cimbala

% Given: Two data sets, A and B of the same variable.
% Goal: Determine whether there is a statistically significant improvement in the variable.

function Two_samples_t_test_same_n
    clear all
    close all

    % First, input data
    
    %{
    % Method 1: Explicit data values
    x_A = [25.6,27.3,24.2,28.7,23.6,25.1]';
    x_B = [26.2,27.1,24.1,29.2,24.5,24.9]';
    %}
    % ---OR---
    
    %{
    % Method 2: Import from Excel
    % NOTE: You must save the Excel file in the same folder as the Matlab m-file
    x_A = xlsread('Two_samples_t_test_same_n.xls',1,'C8:C13');
    x_B = xlsread('Two_samples_t_test_same_n.xls',1,'D8:D13');
    %}
    % ---OR---
    
    % Method 3: Import from text file
    % NOTE: You must save the text file in the same folder as the Matlab m-file
    % Matlab will search the working directory for a file with this name.
    fid1=fopen('Two_samples_t_test_same_n.txt');
    num_lines = 6;
    read_data = textscan(fid1,'%f%f',num_lines);    
    data = [read_data{1,:}];
    x_A = data(:,1);
    x_B = data(:,2);
    
    delta = x_B - x_A;

    avg_A = mean(x_A);
    avg_B = mean(x_B);
    avg_delta = mean(delta);
    
    S_A = std(x_A);
    S_B = std(x_B);
    S_delta = std(delta);
    
    n = length(delta);
    df = n - 1;
    
    fprintf('\n----------------------\n');
    fprintf('Data #\t\tx_A\t\tx_B\t\tdelta = x_B - x_A\n');
    for i = 1:n
        fprintf('%6g\t%7g\t%7g\t\t%8g\n',[i;x_A(i);x_B(i);delta(i)]);
    end
    fprintf('Mean:\t%7g\t%7g\t\t%8g\n',avg_A,avg_B,avg_delta);
    fprintf('S:\t\t%7g\t%7g\t\t%8g\n',S_A,S_B,S_delta);
    fprintf('n = %g\n',n);
    fprintf('df = %g\n',df);
    
    % Null hypothesis: Sample B has a smaller mean than Sample A
    %   OR, mathematically, population mean of variable delta < 0
    
    mu_0 = input('\nEnter mu_0: '); % User input: mu_0
    if isempty(mu_0)
        mu_0 = 0; % Default value
    end
    
    % Calculate the t-statistic and the p-value:
    t = (avg_delta-mu_0)/(S_delta/sqrt(n)); % t-statistic

    tails = 0;
    while (tails ~= 1 && tails ~= 2)
        tails = input('1 or 2 tails? '); % User input: 1-tailed or 2-tailed
        if isempty(tails)
            tails = 1; % Default value
        end
        if (tails ~=1 && tails ~=2) % Error check
            fprintf('Error: You can only use a 1-tailed or 2-tailed test.\n');
        end
    end

    if tails == 2
        fprintf('\nNull hypothesis: mu = %g\n',mu_0);
    else if avg_delta < mu_0
            fprintf('\nNull hypothesis: mu > %g\n',mu_0);
        else if avg_delta > mu_0
                fprintf('\nNull hypothesis: mu < %g\n',mu_0);
            end
        end
    end
    
    p = tdist(t,df,tails); % p-value
    
    fprintf('\nt = %g\n',t);
    fprintf('number of tails = %g\n',tails);
    fprintf('p = %g\n',p);
    
    % Conclusions:
    fprintf('\nConclusion:\n');
    fprintf('There is a %g%% probability that the null hypothesis IS correct.\n',100*p);
    fprintf('There is a %g%% probability that the null hypothesis IS NOT correct.\n',100*(1-p));
    if mu_0 == 0 && avg_delta > 0
        fprintf('We are confident to %g%% that is an increase in the variable.\n',...
            round(100*(1-p)*100)/100);
    else if mu_0 == 0 && avg_delta < 0
            fprintf('We are confident to %g%% that there is a decrease in the variable.\n',...
                round(100*(1-p)*100)/100);
        end
    end
    
    if p <= 0.05
        fprintf('\nSince the confidence level is > 95%%, there is a\n');
        if avg_delta > 0
            fprintf('statistically significant improvement (increase in the mean)\n');
        else if avg_delta < 0
                fprintf('statistically significant change (decrease in the mean)\n');
            end
        end
        fprintf('from Sample A to Sample B.\n');
    else if p >= 0.95
            fprintf('\nTo 95%% confidence, we accept the null hypothesis.\n');
        else
            fprintf('\nSince the confidence level is < 95%%, we cannot claim\n');
            if avg_delta > 0
                fprintf('a statistically significant improvement (increase in the mean)\n');
            else if avg_delta < 0
                    fprintf('a statistically significant change (decrease in the mean)\n');
                end
            end
        end
    end
    site = '"http://www.mne.psu.edu/me345/Lectures/Two_samples_t_test_same_n.xls"';
    title = 'Excel file';
    fprintf(['\nView corresponding <a href = ' site '>' title '</a>\n']);
end

%--------------------------------------------------------------------------

% Function to calculate probability
% Given t-value, degrees of freedom, and number of tails
function p = tdist(t,df,tails)
    if (t < 0)
        t = -t; % Use only positive t
    end
    y = gamma((df+1)/2)/gamma(df/2)/sqrt(pi*df)/(1+t^2/df)^((df+1)/2);
    p = 0.5*betainc(df/(df+t^2),df/2,0.5);
    if (tails == 2)
        p = 2*p; % 2-tailed test option (double probability)
    end
end
