function [Ty,diff,se,CI,p,tcrit,df]=limo_yuen_ttest(a,b,percent,alpha)

% Computes Ty (Yuen's T statistic) to compare the trimmed means of two
% INDEPENDENT groups (single-trials, or subjects with the same number of electrodes).
% Ty=(tma-tmb) / sqrt(da+db), where tma & tmb are trimmed means of a & b,
% and da & db are yuen's estimate of the standard errors of tma & tmb.
% Ty is distributed approximately as Student's t with estimated degrees of freedom, df.
% The p-value (p) is the probability of obtaining a t-value whose absolute value
% is greater than Ty if the null hypothesis (i.e., tma-tmb = mu) is true.
% In other words, p is the p-value for a two-tailed test of H0: tma-tmb=0;
%
% FORMAT [Ty,diff,se,CI,p,tcrit,df]=limo_yuen_ttest(a,b,percent,alpha)
%
% INPUT: a & b are matrices EEGLAB format with trials as the third dimension;
%        percent must be a number between 0 & 100.
%
%   Default values:
%   percent = 20;
%   alpha = 0.05;
%
% OUTPUTS:
%
% Ty = Yuen T statistics
% diff = difference between trimmed means of a and b.
% se standard error
% CI = confidence interval around the difference.
% p = p value.
% tcrit = 1-alpha/2 quantile of the Student's t distribution with adjusted
% df  degrees of freedom
%
% See Wilcox (2005), Introduction to Robust Estimation and Hypothesis
% Testing (2nd Edition), page 159-161 for a description of the Yuen
% procedure.
% You can also check David C. Howell, Statistical Methods for Psychology,
% sixth edition, p.43, 323, 362-363.
%
% See also LIMO_TTEST LIMO_YUEND_TTEST
% Original code by Prof. Patrick J. Bennett, McMaster University
% Various editing, GAR, University of Glasgow, Dec 2007
% Vector version for EEG: GAR, University of Glasgow, Nov 2009
% CI output: GAR, University of Glasgow, June 2010
% ------------------------------
%  Copyright (C) LIMO Team 2019


%% data checking

if nargin<4;alpha=.05;end
if nargin<3;percent=20;end

if size(a,1) ~= size(b,1)
    error('limo_yuen_ttest:InvalidInput', 'Dimension 1 (electrodes) must be the same for a & b');
end

if isempty(a) || isempty(b)
    error('limo_yuen_ttest:InvalidInput', 'DATA cannot have length=0');
end

if (percent >= 100) || (percent < 0)
    error('limo_yuen_ttest:InvalidPercent', 'PERCENT must be between 0 and 100.');
end

if percent == 50
    error('limo_yuen_ttest:InvalidPercent', 'PERCENT cannot be 50, use a method for medians instead.');
end


%% do the computation

if isnan(a(:))
    option = 2;
elseif isnan(b(:))
    option = 2;
else
    option = 1;
end

switch option

    case {1} % there are no nan, no need to loop electrodes
        % -----------------------------------------------------

        % number of trials
        na=size(a,3);
        nb=size(b,3);

        % number of items to winsorize and trim
        ga=floor((percent/100)*na);
        gb=floor((percent/100)*nb);

        % winsorise a
        asort=sort(a,3);
        wa=asort;
        wa(:,:,1:ga+1)=repmat(asort(:,:,ga+1),[1 1 ga+1]);
        wa(:,:,na-ga:end)=repmat(asort(:,:,na-ga),[1 1 ga+1]);
        wva=var(wa,0,3);

        % winsorise b
        bsort=sort(b,3);
        wb=bsort;
        wb(:,:,1:gb+1)=repmat(bsort(:,:,gb+1),[1 1 gb+1]);
        wb(:,:,nb-gb:end)=repmat(bsort(:,:,nb-gb),[1 1 gb+1]);
        wvb=var(wb,0,3);

        % yuen's estimate of standard errors for a and b
        ha=na-2*ga;
        da=((na-1)*wva)/(ha*(ha-1));

        hb=nb-2*gb;
        db=((nb-1)*wvb)/(hb*(hb-1));
        se = sqrt(da+db);
        
        % trimmed means
        ma=mean(asort(:,:,(ga+1):(na-ga)),3);
        mb=mean(bsort(:,:,(gb+1):(nb-gb)),3);

        diff=ma-mb; 

        Ty = diff ./ sqrt(da+db);

        df= ((da+db).^2) ./ ( (da.^2/(ha-1)) + (db.^2/(hb-1)) );

        p=2*(1-tcdf(abs(Ty),df)); % 2-tailed probability

        tcrit=tinv(1-alpha./2,df); % 1-alpha/2 quantile of Student's distribution with df degrees of freedom

        CI(:,:,1) = diff - tcrit.*sqrt(da+db);
        CI(:,:,2) = diff + tcrit.*sqrt(da+db);
        

    case {2} % there are NaNs, so need to loop electrodes
    % ---------------------------------------------------

    tmpa = a;
    tmpb = b;
    
    ma=NaN(size(a,1),size(a,2));
    mb=ma;diff=ma;Ty=ma;df=ma;p=ma;tcrit=ma;
    CI=zeros(size(a,1),size(a,2),2);
    
    for electrode = 1:size(a,1)

        a = squeeze(tmpa(electrode,:,:)); a = a(~isnan(a));
        b = squeeze(tmpb(electrode,:,:)); b = b(~isnan(b));
        
        % number of trials
        na=size(a,2);
        nb=size(b,2);

        % number of items to winsorize and trim
        ga=floor((percent/100)*na);
        gb=floor((percent/100)*nb);

        % winsorise a
        asort=sort(a,2);
        wa=asort;
        wa(:,1:ga+1)=repmat(asort(:,ga+1),[1 ga+1]);
        wa(:,na-ga:end)=repmat(asort(:,na-ga),[1 ga+1]);
        wva=var(wa,0,2);

        % winsorise b
        bsort=sort(b,3);
        wb=bsort;
        wb(:,1:gb+1)=repmat(bsort(:,gb+1),[1 gb+1]);
        wb(:,nb-gb:end)=repmat(bsort(:,nb-gb),[1 gb+1]);
        wvb=var(wb,0,2);

        % yuen's estimate of standard errors for a and b
        ha=na-2*ga;
        da=((na-1)*wva)/(ha*(ha-1));

        hb=nb-2*gb;
        db=((nb-1)*wvb)/(hb*(hb-1));
        se(electrode) = sqrt(da+db);
        
        % trimmed means
        ma(electrode,:)=mean(asort(:,(ga+1):(na-ga)),2);
        mb(electrode,:)=mean(bsort(:,(gb+1):(nb-gb)),2);

        diff(electrode,:)=ma-mb;

        Ty(electrode,:) = diff(electrode,:) ./ sqrt(da+db);

        df(electrode,:) = ((da+db).^2) ./ ( (da.^2/(ha-1)) + (db.^2/(hb-1)) );

        p(electrode,:)=2*(1-tcdf(abs(Ty(electrode,:)),df)); % 2-tailed probability

        tcrit(electrode,:)=tinv(1-alpha./2,df); % 1-alpha/2 quantile of Student's distribution with df degrees of freedom

        CI(electrode,:,1) = diff(electrode,:) - tcrit.*sqrt(da+db);
        CI(electrode,:,2) = diff(electrode,:) + tcrit.*sqrt(da+db);
        
    end

end

end

