% skeweffects.m, based onunequalvariances.m
% Test of effects of positive skew on 1-sample tests.
% Written by Eugene.Gallagher@umb.edu 2001, Revised 3/1/11
% Vary the following parameters
Trials=1000;     % This is a slow routine 400 takes 30 seconds (try 400) 
Xn=100;         % Sample size for X
Xmean=0;    % Mean of X
StdX=1;         % Standard Deviation for X;
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=1;       % To match Sleuth Display 3.4 use skew 1, 0.75 and 0.5
kurtX = 7;     % the normal distribution has kurtosis = 3
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;
TallySign=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)
elseif normal ==0
        fprintf('Uniform random numbers:\n')
        fprintf('Sample X: n=%d  sigma=%3.1f\n',...
            Xn,StdX)
else
        fprintf('Lognormal random numbers:\n')
        fprintf('Sample X: n=%d  sigma=%3.1f\n',...
            Xn,StdX)
end
for i=1:Trials
     if normal==1
        X=pearsrnd(0,StdX,skewX,kurtX,Xn,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
                                % Variance=100^2 * 1/12, see theorem 3.13.1
     else
         X=lognrnd(0,StdX,Xn,1); % Lognormal distribution
     end
    if i==1
        hist(X);
        figure(gcf);pause
    end
    % Perform Wilcoxon's signed rank sum test for difference in location from 0.
%	[pvalueWRS]=wilcoxranksum(X,Y); %Wilcoxon rank sum test, large sample Normal approx
%   pvalueWRS = ranksum(X,Y); % will use exact test for small samples.
    pvalueWSR = signrank(X,0,'alpha',0.05,'method','approximate');
	if pvalueWSR<=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]=ttest(X,0);
	if pvalueSt<=0.05
		TallySt=TallySt+1;
	end
     [pvaluesign,H] = signtest(X,0,'alpha',0.05,'method','exact');
	if pvaluesign<=0.05
		TallySign=TallySign+1;
	end 
end
fprintf(...
 '\nX has %4.0f cases and std(X)=%6.4f\n',...
    length(X),StdX)
fprintf(...
    'There were %d trials\n',...
    Trials);
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)
TallySign=max(1,TallySign); 
fprintf(...
    'Sign test: P(Type I error|alpha=5%%)=%6.3g%%\n',...
    100*TallySign/Trials)