function [results]=tost(sample1,sample2,dL,dU)
% This function performs equivalence testing for 2 independent samples 
% using the TOST procedure (two one-sided tests).
% If BOTH tests are significant (p<0.05), the two samples are equivalent.
% In this case, report the statistics for the test with the HIGHEST p-value. 
%-----------------------------------------------------------------------------
% Inputs: sample 1, sample 2, lower equivalence bound in Cohen's d, upper
% equivalence bound in Cohen's d.
%-----------------------------------------------------------------------------
% Output: "results" structure with the following fields:
% 1) t-value and p-value for the lower equivalence bound (one-sided t-test #1)
% 2) t-value and p-value for the upper equivalence bound (one-sided t-test #2)
% 3) degrees of freedom
%-----------------------------------------------------------------------------
% Reference: 
% Lakens, D. (2017). Equivalence tests: a practical primer for t tests, 
% correlations, and meta-analyses. Social Psychological and Personality 
% Science, 8(4), 355-362.
%-----------------------------------------------------------------------------

n1=sum(~isnan(sample1));
n2=sum(~isnan(sample2));
M1=nanmean(sample1);
M2=nanmean(sample2);
SD1=nanstd(sample1);
SD2=nanstd(sample2);

sd_pooled=sqrt(((n1-1)*SD1^2+(n2-1)*SD2^2)/(n1+n2-2));

% Convert Cohen's d bounds to raw scores.
dL=dL*sd_pooled;
dU=dU*sd_pooled;

% Calculate the t-values and df's
if n1==n2 % Student's t-test when the number of observations in the 2 samples is the same
    tL=((M1-M2)-dL)/(sd_pooled*sqrt((1/n1)+(1/n2)));
    tU=((M1-M2)-dU)/(sd_pooled*sqrt((1/n1)+(1/n2)));
    df=n2+n1-2;
else % Welch's t-test when the number of observations in the 2 samples is different
    tL=((M1-M2)-dL)/(sqrt((SD1^2/n1)+(SD2^2/n2)));
    tU=((M1-M2)-dU)/(sqrt((SD1^2/n1)+(SD2^2/n2)));
    numer=((SD1^2/n1)+(SD2^2/n2))^2;
    denom=((SD1^2/n1)^2/(n1-1))+((SD2^2/n2)^2/(n2-1));
    df=numer/denom;
end

% Calculate the p-values
if tL>0 
    pL=1-tcdf(abs(tL),df);
else 
    pL=tcdf(abs(tL),df);
end
if tU<0 
    pU=1-tcdf(abs(tU),df);
else 
    pU=tcdf(abs(tU),df);
end

results.tL=tL;
results.tU=tU;
results.df=df;
results.pL=pL;
results.pU=pU;

end

