% unequalvariances.m
% Test of effects of unequal variances on 2-sample tests.
% Written by Eugene.Gallagher@umb.edu 2001, Revised 11/22/2010,2/02/11
% Vary the following parameters
Trials=200;     % This is a slow routine 400 takes 30 seconds (try 400)
Trialsrp=100;   % Trials for each random permutation test (try 400) 
Xn=4;         % Sample size for X
Xmean=0;    % Mean of X
StdX=1;         % Standard Deviation for X;
Yn=20;          % Sample size for Y
Ymean=0;      % Mean of Y
StdY=.5235;         % Standard Deviation for Y;
normal=1;       % 1 for normal; Change to 0 for a uniform distribution; 
%                 2 for lognormal
% skew & kurt only defined if normal distribution is specified.
skewX=0;          % normal distribution has skewness, 0
skewY=0;          % To match Sleuth Display 3.4 use skew 1, 0.75 and 0.5
kurtX = 3;        % the normal distribution has kurtosis = 3
kurtY = 3;         % To match sleuth use kurtosis 6 and kurtosis 2.5
UseT=1;         % 1 if t test, 0 if |difference| for random permutations
% Create some variables to hold the results of the test
% if the p value is less than 0.05, then the tally will be increased by 1
% P(Type I error), or alpha is set at 0.05 for the tests, so Tally*/Trials
% should be 5% if the test is working properly.
TallySt=0;
TallyWelch=0;
TallyWilcox=0;
TallyWilcoxE=0;
Tallyrandp=0;
     % print out the parameters:
if normal ==1
        fprintf('Pearson normal random numbers:\n')
        fprintf('Sample X: n=%d Mean=%7.3g sigma=%3.1f, skew=%4.2f, kurt=%4.2f\n',...
            Xn, Xmean,StdX,skewX,kurtX);
        fprintf('Sample Y: n=%d Mean=%7.3g sigma=%3.1f, skew=%4.2f, kurt=%4.2f\n',... 
            Yn,Ymean,StdY,skewY,kurtY);
elseif normal ==0
        fprintf('Uniform random numbers:\n')
        fprintf('Sample X: n=%d  sigma=%3.1f\n',...
            Xn,StdX)
        fprintf('Sample Y: n=%d  sigma=%3.1f\n',... 
            Yn,StdY)
else
        fprintf('Lognormal random numbers:\n')
        fprintf('Sample X: n=%d  sigma=%3.1f\n',...
            Xn,StdX)
        fprintf('Sample Y: n=%d  sigma=%3.1f\n',... 
            Yn,StdY)
end
for i=1:Trials
     if normal==1
        X=pearsrnd(0,StdX,skewX,kurtX,Xn,1);
        Y=pearsrnd(0,StdY,skewY,kurtY,Yn,1);
     elseif normal==0
         X=StdX*sqrt(12)*(rand(Xn,1)-0.5);% a uniform distribution variance 1/12th,
				                % mean 0., see LMex040302.m for the 1/12
         Y=StdY*sqrt(12)*(rand(Yn,1)-0.5);  % a uniform distribution, 
                                % Variance=100^2 * 1/12, see theorem 3.13.1
     else
         X=lognrnd(0,StdX,Xn,1); % Lognormal distribution
         Y=lognrnd(0,StdY,Yn,1);
     end
    if i==1
        subplot(2,1,1)
        hist(X);
        subplot(2,1,2)
        hist(Y);
        figure(gcf);pause
        subplot(1,1,1)
    end
    % Perform Wilcoxon's rank sum test for difference in median between X & Y.
%	[pvalueWRS]=wilcoxranksum(X,Y); %Wilcoxon rank sum test, large sample Normal approx
%   pvalueWRS = ranksum(X,Y); % will use exact test for small samples.
    pvalueWRS = ranksum(X,Y,'alpha',0.05,'method','approximate');
	if pvalueWRS<=0.05
		TallyWilcox=TallyWilcox+1;
	end
    % Perform Student's t test for difference in mean. The means are set to be the same.
%	[D,t,df,pvalueSt,CI]=stud2sample(X,Y);
     [H,pvalueSt]=ttest2(X,Y);
	if pvalueSt<=0.05
		TallySt=TallySt+1;
	end
	% the unequal variance t test, also called the Welch test
	% [D,t,df,pvalueWelch,CI]=welch2sample(X,Y);
     [H,pvalueWelch] = ttest2(X,Y,0.05,'both','unequal');
	if pvalueWelch<=0.05
		TallyWelch=TallyWelch+1;
	end 
	% This is the 2-sample random permutation test, which also
	% also may be affected by unequal variances.
    % This random permutation approach slows down the trials incredibly.
    % Undelete if you have a fast processor
    %  [D,Studp,t,df,RandPermP]=randp2sample(X,Y,Trialsrp);
    if UseT==1
         RandPermP=randpermtratio(X,Y,Trialsrp);
    else
        [prob,meandiff,tally]=randpermmean(X,Y,Trialsrp);
    end
	if RandPermP<=0.05
		 Tallyrandp=Tallyrandp+1;
    end
end
fprintf(...
 '\nX has %4.0f cases; Y has %4.0f cases;\nstd(X)=%6.4f; std(Y)=%6.4f\n',...
    length(X),length(Y),StdX,StdY)
fprintf(...
    'There were %d trials with %d random permutations within each trial.\n',...
    Trials,Trialsrp);
TallySt=max(1,TallySt);
fprintf('Student''s t: P(Type I error| alpha=5%%)=%6.3g%%\n',...
    100*TallySt/Trials);
TallyWilcox=max(1,TallyWilcox);
fprintf(...
'Wilcoxon''s rank sum test, approx: P(Type I error|alpha=5%%)=%6.3g%%\n',...
    100*TallyWilcox/Trials)
TallyWelch=max(1,TallyWelch); 
fprintf(...
    'Welch''s unequal variance t test: P(Type I error|alpha=5%%)=%6.3g%%\n',...
    100*TallyWelch/Trials)
TallySt=max(1,TallySt);
if UseT==1
    fprintf('Random permutations of t ratios: P(Type I error|alpha=5%%)=%6.3g%%\n',...
    100*TallySt/Trials)
else
    fprintf('Random permutations of means: P(Type I error|alpha=5%%)=%6.3g%%\n',...
    100*TallySt/Trials)
end