function [ci,tstat,nu,pvalue]=f(x,y, percent,delta)
% twosampleci(x,y,percent,delta) computes Welch's approximate confidence interval for the difference of
% of the means mu_x  - mu_y (which can also be used for a two-sided test)
% at the confidence "percent" in the range 0-100.
% x, y are two vectors of observations, possibly of different lengths.
% It also gives the value tstat of the t-statistic for testing mu_x-mu_y=delta and the two-sided p-value.
% For paired data, use the one-sample method.
if nargin < 3, 
percent=95; 
end
if nargin < 4, delta=0; 
end
sx=std(x);
sy=std(y);
m=length(x);
n=length(y);
nu=(sx^2/m + sy^2/n)^2/( (sx^2/m)^2/(m-1) + (sy^2/n)^2/(n-1) );
tc=tinv(1- (100-percent)/200,nu); 
l=mean(x)-mean(y) - tc*sqrt(sx^2/m + sy^2/n);
u=mean(x)-mean(y) + tc*sqrt(sx^2/m + sy^2/n);
ci=[l; u];
tstat=(mean(x)-mean(y) -delta)/sqrt(sx^2/m + sy^2/n);
pvalue=2*tcdf(-abs(tstat),nu);