Not working
function [D,pvalue,CI]=MC2sample(X,Y,alpha,Trials)
% [D,pvalue,CI]=MC2sample(X,Y,alpha,Trials)
% Input: X,Y two samples
%        alpha level for CI, optional, 95% if not specified
%        Number of Monte Carlo trials
% Output: D=difference in mean;
%         p value, 2-sided
%         CI for 1-alpha CI, [L U];95% CI if alpha not specified.
% Written by E. Gallagher;
X=X(:);Y=Y(:);
D=mean(X)-mean(Y);
n=length(X);
m=length(Y);
XY=[X;Y];
sumXY=sum(XY(:,1));
if nargin<4
    Trials=1000;
end
Tally=zeros(Trials,1);
for i=1:Trials
    j=randperm(n+m);
    Tally(i)=mean(XY(j(1:n)))-mean(XY(j(n+1:n+m)));
end
i=length(find(abs(Tally)>=abs(D)));
pvalue=i/Trials;

if nargin<4;
% Use stixbox's pt.m for significance of t.
   pvalue=2*(1-pt(t,df));
end
if nargin<3;
    alpha=0.05;
end
% Use qt.m, the inverse Student's t distribution
% Theorem 9.5.1 for confidence interval.
HalfCI=qt(1-alpha/2,df)*Sp*sqrt(1/n+1/m);
CI=[D-HalfCI D+HalfCI];