% Two Samples - Matlab t-test for Two Independent Samples
% Case with different values of n
% 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_different_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,26.5,25.8]';
    %}
    % ---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_different_n.xls',1,'C8:C13');
    x_B = xlsread('Two_samples_t_test_different_n.xls',1,'D8:D15');
    %}
    % ---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_different_n.txt');
    num_lines = 8;
    fseek(fid1,0,'bof');
    read_data = textscan(fid1,'%f%f',num_lines);    
    data = [read_data{1,:}];
    x_A = data(:,2);
    x_B = data(:,1);
    x_A = x_A(1:end-sum(isnan(x_A))); % remove NaN from data
    x_B = x_B(1:end-sum(isnan(x_B)));
    
    n_A = length(x_A); % Sample number of data points (n)
    df_A = n_A - 1; % Degrees of freedom (df)
    avg_A = mean(x_A); % Sample mean (x_bar)
    S_A = std(x_A); % Sample standard deviations (S)
    var_A = S_A^2; % Sample variances (S^2)
    
    n_B = length(x_B); % Sample number of data points (n)
    df_B = n_B - 1; % Degrees of freedom (df)
    avg_B = mean(x_B); % Sample mean (x_bar)
    S_B = std(x_B); % Sample standard deviations (S)
    var_B = S_B^2; % Sample variances (S^2)
    
    % We use Welch's equation to calculate the degrees of freedom:
    df_real = ((var_A/n_A+var_B/n_B)^2)/((var_A/n_A)^2/df_A+(var_B/n_B)^2/df_B); 
    df = round(df_real); % degrees of freedom (nearest integer)
    
    fprintf('\n----------------------\n');
    fprintf('Data #\t\tx_A\t\tx_B\n');
    for i = 1:min(n_A,n_B)
        fprintf('%6g\t%7g\t%7g\n',[i;x_A(i);x_B(i)]);
    end
    for i = min(n_A,n_B)+1:max(n_A,n_B)
        if n_A == min(n_A,n_B)
            fprintf('%6g\t\t\t%7g\n',[i;x_B(i)]);
        else
            fprintf('%6g\t%7g\n',[i;x_A(i)]);
        end
    end
    fprintf('n = \t%7g\t%7g\n',n_A,n_B);
    fprintf('df = \t%7g\t%7g\n',df_A,df_B);
    fprintf('Mean:\t%7g\t%7g\n',avg_A,avg_B);
    fprintf('S:\t\t%7g\t%7g\n',S_A,S_B);
    
    % Null hypothesis: Sample B has a smaller mean than Sample A
    %   OR, mathematically, mu_B < mu_A
    
    % Calculate the t-statistic and the p-value:
    t = (avg_A-avg_B)/sqrt(var_A/n_A+var_B/n_B); % 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_B = %g\n',avg_A);
    else if avg_B < avg_A
            fprintf('\nNull hypothesis: mu_B > %g\n',avg_A);
        else if avg_B > avg_A
                fprintf('\nNull hypothesis: mu_B < %g\n',avg_A);
            end
        end
    end
    
    p = tdist(t,df,tails); % p-value
    
    fprintf('\nt = %g\n',t);
    fprintf('df_real = %g\n',df_real);
    fprintf('df = %g\n',df);
    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 tails == 2
        fprintf('We are confident to %g%% that there is no change in the variable.\n',...
            round(100*p*100)/100);
    else if avg_B > avg_A
        fprintf('We are confident to %g%% that is an increase in the variable.\n',...
            round(100*(1-p)*100)/100);
        else if avg_B < avg_A
            fprintf('We are confident to %g%% that there is a decrease in the variable.\n',...
                round(100*(1-p)*100)/100);
            end
        end
    end
    
    if p <= 0.05
        fprintf('\nSince the confidence level is > 95%%, there is a\n');
        if tails == 2
            fprintf('statistically significant result that there is no difference in the means\n');
        else if avg_B > avg_A
                fprintf('statistically significant improvement (increase in the mean)\n');
                else if avg_B < avg_A
                        fprintf('statistically significant change (decrease in the mean)\n');
                    end
            end
            fprintf('from Sample A to Sample B.\n');
        end
    else if p >= 0.95
            fprintf('\nTo 95%% confidence, we accept the null hypothesis that mu_A = mu_B.\n');
        else
            fprintf('\nSince the confidence level is < 95%%, we cannot claim\n');
            if tails == 2
                fprintf('that the two means are equal.\n');
            else if avg_B > avg_A
                    fprintf('a statistically significant improvement (increase in the mean).\n');
                    else if avg_B < avg_A
                            fprintf('a statistically significant change (decrease in the mean).\n');
                        end
                end
                fprintf('from Sample A to Sample B.\n');
            end
        end
    end
    site = '"http://www.mne.psu.edu/me345/Lectures/Two_samples_t_test_different_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
