function [D,t,df,pvalue,CI]=Welch2sample(X,Y,alpha)
% Student't 2-sample unequal variance t test.
% [D,t,df,pvalue,CI]=Welch2sample(X,Y,alpha)
% Input: X,Y two samples
%        alpha level for CI, optional, 95% if not specified
% Output: D=mean(X)-mean(Y);
%         t=Student's t statistic
%         df=degrees of freedom for t statistic
%         p value, 2-sided, for t Statisitc, with df degrees of freedom
%         CI for 1-alpha CI, [L U];95% CI if alpha not specified.
% Written by E. Gallagher; Revised 11/14/10
X=X(:);Y=Y(:);
n=length(X);
m=length(Y);
D=mean(X)-mean(Y);
% Estimated standard error (using separate estimates of variance for each sample)

SEW=sqrt(var(X)/n+var(Y)/m);  %Sleuth p. 97
t=D/SEW;
% Saterthwaite's approximation for d.f.
df=floor(SEW.^4/((std(X)/sqrt(n)).^4 /(n-1) + (std(Y)/sqrt(m)).^4 /(m-1)));
% Use tcdf from Matlab's statistics toolbox:
if t>=0
    pvalue = 2*(1-tcdf(t,df));
else
    pvalue=2*tcdf(t,df);
end
if nargin<3;
    alpha=0.05;
end

HalfCI=tinv(1-alpha/2,df)*SEW;
CI=[D-HalfCI D+HalfCI];