%%%this class has output= HotelingtTest(self,M, withinOrBetween), in output
%%% output=[tsquared, df1_tsquared, df2_tsquared, F, df1,df2, p, d]; note
%%% that d is dz for within and d for between
%This class also has [n, df2, critical_F, F, lambdaF, p2]=Hotelingtpowerd(self,d,k,n,alpha,withinOrBetween)
% note that d is define as t/sqrt(2/n), so if we use within and dz, use
% d=dz*sqrt(2)
%this class has lr01= HotelingtlrfromDs(self,dob,dtie,withinOrBetween,k,n)
%it is based on an point alternative, dob and dtie are both defined by t/sqrt(2/n)
%BF01=HotelingBF01BasedonD_cauchy(self,dob, k, n, withinOrBetween,scale)
%it is based on a cauchy alternative, dob is defined by t/sqrt(2/n)
% note that all these functions can be used for any dimension data, k=1 is
% univariate
classdef collectiveHotelingTsquared
    methods
        function [output,raw]= HotelingtTest(self,M, withinOrBetween)
        %input parameters: M should be nX(2k) matrix, the second half collumns would be the same values across rows if it is used for one sample test, withinOrBetween=1 for within (for paired or one sample), 2 for between (indepdent sample)
        %output tsquared is hoteling tsquared, F is, p is p value, 
        % d is effect size, tsquared=d^2*n/2 when withinOrBetween=2 (between), tsquared=d^2*n when withinOrBetween=1 (within, like dz),
            k=length(M(1,:))/2; %k is the dimension of multivariate. k=1 is univariate
            n=length(M(:,1));
            if withinOrBetween==1
                Md=[]
                for i=1:k
                    Md=[Md,M(:,i)-M(:,i+k)];
                end 
        
                %Md=[M(:,1)-M(:,3),M(:,2)-M(:,4)];
                COVMd=cov(Md);%unbiased cov
                MeanMd=mean(Md);
                
                dsquared=MeanMd*inv(COVMd)*MeanMd';
                d=sqrt(dsquared); %d is similar to cohen's dz
                tsquared=d^2*n;     

                rawd=MeanMd;
                rawcov=COVMd;
            elseif withinOrBetween==2
                %M1=[M(:,1),M(:,2)];
                %M2=[M(:,3),M(:,4)];
                M1=M(:,1:k);
                M2=M(:,k+1:end);
                MeanM1=mean(M1);
                MeanM2=mean(M2);
                COVM1=cov(M1);
                COVM2=cov(M2);
                COVpool=(COVM1+COVM2)/2;
                %n=length(M(:,1));
                dsquared=(MeanM1-MeanM2)*inv(COVpool)*(MeanM1-MeanM2)';
                d=sqrt(dsquared); %d is similar to cohen's d
                tsquared=d^2*n/2;
                rawd=MeanM1-MeanM2;
                rawcov=COVpool;
            else
                print('withinOrBetween should be 1 or 2')
            end
            
            df1_tsquared=k;
            df2_tsquared=withinOrBetween*(n-1);
        
            F=tsquared*(df2_tsquared-df1_tsquared+1)/(df1_tsquared*df2_tsquared); % F(2,n-2)
            df1=df1_tsquared;
            df2=df2_tsquared-df1_tsquared+1;
            p=1-ncfcdf(F,df1,df2,0);
            % output.Hotelingtsquare=tsquared;
            % output.d=d;
            % output.F=F;
            % output.p=p;

            

            output=[tsquared, df1_tsquared, df2_tsquared, F, df1,df2, p, d];
            raw.d=rawd,
            raw.cov=rawcov;
        end
        function  [n, df2, critical_F, F, lambdaF, p2]=Hotelingtpowerd(self,d,k,n,alpha,withinOrBetween)
        %d=t*sqrt(2/N), which is sqrt(2)*dz for within design; withinOrBetween=1 for within, 2 for between
        %k is the diension of data 
            % here d is defined as following
            % Hoteling tsquared
            % dsquared=d^2
            % tsquared = dsquared*n/2
            % so dz=dsquared/sqrt(2) if we use dz in within situation and d in
            % between case
            
            tsquared=d^2*n/2;
            df1_tsquared=k;
            df2_tsqaured=withinOrBetween*(n-1);
        
            F=tsquared.*(df2_tsqaured-df1_tsquared+1)./(df1_tsquared.*df2_tsqaured); % F(2,n-2)
            df1=df1_tsquared;
            df2=df2_tsqaured-df1_tsquared+1;
            lambdaF=tsquared; %lambdaF is sqaure of predicted t or predicted hoteling Tsquared
            %lambdaF1=F*df1;
            %critical_F=finv(0.95,df1,df2);
            critical_F=finv(1-alpha,df1,df2);
        
        
            p2=1-ncfcdf(critical_F,df1,df2,lambdaF);
        
            %[n, df1,df2, critical_F, lambdaF, p2]
            
            showFigure=0;
            if showFigure==1
                figure 
                plot(n,p2,'-')
            end
        end
        function lr01= HotelingtlrfromDs(self,dob,dtie,withinOrBetween,k,n)
            % alternative is a point
            % lr=lrfromDs(0.5,0.6585, 2, 20)
            %<<lr=0.319
            % replication study n =20, observed d=0.5, theoratical interesting estimate d=0.6585, between design a=2
            %n=20
            %a=2 %2 for between design, 1 for within design
        
            %lrfromDs(0.25,0.5806, 1,50)% n =50, observed d=0.25, dtie=0.5806, between design a=2
            %<<1.7716
            
            tsquaredob=dob^2*n/2;
            tsquaredtie=dtie^2*n/2; % used for the noncentral F 
            lambdaFtie=tsquaredtie;
            df1_tsquared=k;
            df2_tsquared=withinOrBetween*(n-1);
            
            [Fob,df1,df2]=self.tsquaredtoF(tsquaredob,df1_tsquared,df2_tsquared);
            
            L0=ncfpdf(Fob,df1,df2,0);
            L1=ncfpdf(Fob,df1,df2,lambdaFtie);
            lr01=L0./L1;
            
            %F is related to two tail T, so F LR is the average of two LRs based on
            %t test with the same abs(d), below is to that only
            tob=sqrt(tsquaredob);
            tie=sqrt(tsquaredtie);
            lambdatie=tie;
            L0t=nctpdf(tob,df2_tsquared,0);
            L1t1=nctpdf(tob,df2_tsquared,lambdatie);
            L1t2=nctpdf(tob,df2_tsquared,-lambdatie);
        
            lr01t=2*L0t./(L1t1+L1t2);
          
        end
        
        function BF01=HotelingBF01BasedonD_cauchy(self,dob, k, n, withinOrBetween,scale) 
            %observedD is defined as observed diff/MSE or t*sqrt(2/n); n is subject number in each condition; scale is 1/sqrt(2) by default, withinOrBetter=2 for between or 1 for within
            %k is dimension, k=1 is univariate like BF01BasedonD_cauchy(observedD, n, withinOrBetween,scale) 
        
            %between example
            %MSE=3079.755;		
            %n=20;
        
            %observedt=(69.514-64.65)/sqrt(2*MSE/n);
            %observedD=(69.514-64.65)/sqrt(MSE);
        
            %withinOrBetween=2 %2for between 1 for within
        
            %BF01=HotelingBF01BasedonD_cauchy(observedD, 1, n, withinOrBetween,1/sqrt(2)); 
          %BF01 = 3.1407   
            %within example
            
        %     MSE=470.719;		
        %     n=30;
        % 
        %     observedt=(26.569-19.624)/sqrt(2*MSE/n);
        %     observedD=observedt*sqrt(2/n);
        %     withinOrBetween=1 %2for between 1 for within
        % 
        % 
        %     BF01=HotelingBF01BasedonD_cauchy(observedD, 1, n, withinOrBetween,1/sqrt(2)); 
        %BF01 =    2.5658
            
                
                tsquaredob=dob^2*n/2;
               
                df1_tsquared=k;
                df2_tsquared=withinOrBetween*(n-1);
                
                [Fob,df1,df2]=self.tsquaredtoF(tsquaredob,df1_tsquared,df2_tsquared);
                
                pd = makedist('tLocationScale','mu',0,'sigma',scale,'nu',1); %@(x)(1/(pi*scale)*(scale^2./(x.^2+scale^2))) it is cauchy PDF to get the pdf of d for between design and dz for within design
                %fun2=@(x)(1/(pi*scale)*(scale^2./(x.^2+scale^2))).*(nctpdf(observedt,df,x*sqrt(n/2)));
        
                fun2=@(x)pdf(pd,x).*(ncfpdf(Fob,df1,df2,x.^2*n/withinOrBetween)); % Likelihood weighted by pdf of d
                L1 = integral(fun2,-Inf,Inf); %weighted average of the likelihood of the observed Hoteling T squared given alternative Hoteling T squared which is based on the d or dz from cauchy 
                
                L0=ncfpdf(Fob,df1,df2,0);
        
                BF01=L0/L1;
        
        end
        function [F,df1,df2]=tsquaredtoF(self, tsquared,df1_tsquared,df2_tsquared)
            F=tsquared.*(df2_tsquared-df1_tsquared+1)./(df1_tsquared.*df2_tsquared); % F(2,n-2)
            df1=df1_tsquared;
            df2=df2_tsquared-df1_tsquared+1;    
        end
    end
end