function [D,Studp,t,df,RandPermP,CIpD,CIr]=randp2sample(X,Y,Trials,UseT,alpha);
% Test for difference between 2 means using random permutations
% format [D,Studp,t,df,RandPermP,CIpD,CIr]=randp2sample(X,Y,Trials,UseT,alpha);
%           Example: [D,Studp,t,df,RandPermP,CIp,CIr]=randp2sample(X,Y);
% Input:    X, Y  column or row vectors with cases as rows, variables as columns.
%           Trials, number of random permutations [Optional]
%               if nargout>2 & nargin<3, Trials=1e4
%           UseT, ==1, Use t statistic for random permuation tests
%                 ==0, Use abs(D from randperm)>abs(D) for random permutation tests
%                 Note, these two methods produce nearly identical results.
% Output:   D   Difference in means
%           Studp  parametric p value for difference based on Student's t, two-tailed
%           T   Student's t statistic, n+m-2 df
%           df for 2-sample, equal mean t test, equal variance assumed
%           RandPermP   2-tailed p test of D=0, based on Random permutations
%           CIpD Parametric 1-alpha% CI for D based on Student's t
%           CIr 1-alpha% CI for difference of X & Y, 
%           note that this CI is for a difference of 0. It is not the CI for 
%           the observed difference 
% Reference: Legendre & Legendre, Manly
% Written by Eugene.Gallagher@umb.edu
% see also

if nargout>4 & nargin<3
    Trials=1e4;
    UseT=1;
    alpha=0.05;
elseif nargin<5
    alpha=0.05;
end
X=X(:);
n=length(X);
Y=Y(:);
m=length(Y);
[D,t,df,Studp,CIpD]=stud2sample(X,Y,alpha);
df=n+m-2;
XY=[X;Y];
if nargin>3
    UseT=logical(UseT);     % Simple 1's and 0's are evaluated logically, but this is
                            % more formal, and perhaps faster.
else
    UseT=logical(0);
end 
if nargout>3
    Tally=0;   % Create a matrix of all zeros to hold results
    
    if nargout>6
        TallyD=zeros(Trials,1);
    end
    for i=1:Trials
        j=randperm(n+m);
        RandD=mean(XY(j(1:n)))-mean(XY(j(n+1:n+m)));
        if nargout>6
            TallyD(i)=RandD;
        end
        Sp=sqrt(    ((n-1)*var(XY(j(1:n)))+(m-1)*var(XY(j(n+1:n+m))))...
                /  (n+m-2) );
        randt=RandD/(Sp*sqrt(1/n+1/m));
        if UseT     % Use t statistic to determine p value
            if abs(randt)+eps >= abs(t)
                Tally=Tally+1;
            end
        else
            if abs(RandD)+eps>=abs(D) % Two-sided test using absolute value of D
                Tally=Tally+1;
            end
        end
    end
end
RandPermP=Tally/Trials;
if nargout>6
    TallyD=sort(TallyD);
    CIr=zeros(1,2);
    CIr(1)=TallyD(floor(alpha/2*Trials));
    CIr(2)=TallyD(ceil((1-alpha/2)*Trials));
end