% Example - Using the Student's t Distribution to
% Estimate a Mean Value and Confidence Interval

% Barukyah Shaparenko and John Cimbala

function confidence_interval

    clear all
    close all

    % Import data from Excel
    % Save Excel data in same folder as Matlab m-file
    % Use xlsread or xlsxread

    % sheet_number and data_range are optional
    % Default is Sheet 1 with all (numerical) data imported

    % Option 1 - read data from Excel:
    % xls_filename = 'Using_t_PDF_to_estimate_mean.xls';
    % sheet_number = 1;
    % data_range = 'B5:B9';
    % x = xlsread(xls_filename,sheet_number,data_range);
   
    % Option 2 - read data here directly:
    x = [312.2; 320.6; 315.5; 314.1; 319.6];

    % Basic statistics for this data set:

    n = length(x);
    sum_3 = 0;
    sum_4 = 0;
    for i = 1:n
      del_x_avg(i) = x(i) - mean(x);
      sum_3 = sum_3 + del_x_avg(i)^3;
      sum_4 = sum_4 + del_x_avg(i)^4;
    end

    Mean = mean(x);
    Standard_Error = std(x)/sqrt(n);
    Median = median(x);
    Mode = mode(x);
    Standard_Deviation = std(x);
    Sample_Variance = Standard_Deviation^2;
    Kurtosis = n*(n+1)/(n-1)/(n-2)/(n-3)*sum_4/Standard_Deviation^4 - 3*(n-1)^2/(n-2)/(n-3);
    Skewness = n/(n-1)/(n-2)*sum_3/Standard_Deviation^3;
    Range = max(x) - min(x);
    Minimum = min(x);
    Maximum = max(x);
    Sum = sum(x);
    Count = length(x);

    fprintf('\n');
    fprintf('Mean \t\t\t\t%g\n',Mean);
    fprintf('Standard Error \t\t%g\n',Standard_Error);
    fprintf('Median \t\t\t\t%g\n',Median);
    fprintf('Mode \t\t\t\t%g\n',Mode);
    fprintf('Standard Deviation \t%g\n',Standard_Deviation);
    fprintf('Sample Variance \t%g\n',Sample_Variance);
    fprintf('Kurtosis \t\t\t%g\n',Kurtosis);
    fprintf('Skewness \t\t\t%g\n',Skewness);
    fprintf('Range \t\t\t\t%g\n',Range);
    fprintf('Minimum \t\t\t%g\n',Minimum);
    fprintf('Maximum \t\t\t%g\n',Maximum);
    fprintf('Sum \t\t\t\t%g\n',Sum);
    fprintf('Count \t\t\t\t%g\n',Count);

    df = n - 1;
    alpha = 0.05;
    t = tinv(alpha,df);
    S = Standard_Deviation;
    x_bar = Mean;
    mu_min = Mean - t*S/sqrt(n);
    mu_max = Mean + t*S/sqrt(n);
    
    fprintf('\nPart(a)\n');
    fprintf('n = \t\t\t%g\t\tnumber of data points\n',n);
    fprintf('df = \t\t\t%g\t\tdegrees of freedom = n-1\n',df);
    fprintf('alpha = \t\t%g\t(%g percent confidence level)\n',alpha,100*(1-alpha));
    fprintf('t_alpha/2 = \t%g\tUse tinv function - students''t\n',t);
    fprintf('S = \t\t\t%g\tSample standard deviation (from above)\n',S);
    fprintf('x_bar = \t\t%g\tSample mean (from above)\n',x_bar);
    fprintf('mean_min = \t\t%g\tMinimum predicted population mean\n',mu_min);
    fprintf('mean_max = \t\t%g\tMaximum predicted population mean\n',mu_max);
    fprintf('\nInterval: %g < pop. mean < %g\n\n',mu_min,mu_max);
    
    limit = 2;
    
    % Iteration procedure to estimate the required n:
    i = 1;
    n_guess(i) = 10;
    while (i<=10)
        df_guess(i) = n_guess(i) - 1;
        t_guess(i) = tinv(alpha,df_guess(i));
        n_guess(i+1) = (S*t_guess(i)/limit)^2;
        i = i + 1;
    end
    
    fprintf('\nPart(b)\n');
    fprintf('+/- limit = \t%g\tSpecify the confidence interval limit (one side)\n',limit);
    fprintf('\nIteration procedure to estimate the required n:\n');
    fprintf('   Guess n\t\t\tdf=n-1\t\t  t_alpha/2\t\t     new n (n = S^2*t^2/limit^2)\n');
    for k = 1:i-1
        fprintf('%10.4f\t\t%10.4f\t\t%10.4f\t\t%10.4f\n',n_guess(k),df_guess(k),t_guess(k),n_guess(k+1));
    end
    fprintf('\nConverged to n = %g\n',ceil(n_guess(i)));
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

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

% Function to find t
% Given p-value, degrees of freedom
% Uses 2-tailed test automatically
function t = tinv(p,df)
    if (p > 0.5)
        p = 1-p; % Find only positive t
    end
    % Guess initial range
    t_low = 0;
    exponent = -1;
    t_high = inf;
    while(exponent > -10) % Smallest error in t: 10^-10
        t_difference = 10^exponent;
        for t_guess = t_low:t_difference:t_high
            prob = tdist(t_guess,df,2); % Use tdist to find probability
            if prob < p % Compare probability to stated value
                break;
            end
        end
        % Specify new range for t
        t_low = t_guess - t_difference;
        t_high = t_guess;
        exponent = exponent - 1; % Refine search by factor of 10
    end
    t = t_guess;
end
