
function output=ellipsePandE20231125(rawormean, confidencelevel,XY,mark,color) %1 for raw 2 for mean, if confidencelevel > 1 then only draw the 1 std or standard error circle
     %%%%this function draw the residual ellipse with Mu=(0,0),
     %%%%Sigma=[X.^2,X.*Y;X.*Y,Y.^2]

            X=XY;
            Mu = mean( X);
            X0 = bsxfun(@minus, X, Mu);

        %     STD = 2;                     %# 2 standard deviations
        %     conf = 2*normcdf(STD)-1;     %# covers around 95% of population
             conf = confidencelevel; %percentage
             
             if rawormean==1
                 if conf>=1 || conf < 0
                    scale =1;
                 else 
                     scale = chi2inv(conf,2);     %# inverse chi-squared with dof=#dimensions, sqrt of scale is times of standard deviation of x and y
                 end
                 Cov = cov(X0) * scale;
                 %Sigma=[mean(XY(:,1).^2),mean(XY(:,1).*XY(:,2));mean(XY(:,1).*XY(:,2)),mean(XY(:,2).^2)]*scale;
                 
             elseif rawormean==2
                 if conf>=1 || conf < 0
                    scale1 =1;
                 else 
                     scale1 = 2*(size(X,1)-1)/(size(X,1)-2)*finv(conf,2, size(X,1)-2);     %# scale is Hotelling T square, Hotelling T square  = 2(n-1)/(n-2)*F (2,n-2), inverse F (2,n-2) with dof=#dimensions and n-2
                 end
                     %Cov1 = cov(X0) * scale;
                 Cov = cov(X0) * scale1/size(X,1);
                  %Sigma=[mean(XY(:,1).^2),mean(XY(:,1).*XY(:,2));mean(XY(:,1).*XY(:,2)),mean(XY(:,2).^2)]* scale1/size(X,1);
             else   
             end
             
             %Cov=Sigma;
             
             [V, D] = eig(Cov); %# the sum of D is the sum of variances in Cov, which means after sum of variance does not change after eigen decomposition or canonical transformation.
        %     %# eigen decomposition [sorted by eigen values]
        %    [V, D] = eig( X0'*X0 ./ (sum(idx)-1) );     %#' cov(X0)
            [D, order] = sort(diag(D), 'descend');
            D = diag(D);
            V = V(:, order);

            t = linspace(0,2*pi,100);
            e = [cos(t) ; sin(t)];        %# unit circle
             %http://en.wikipedia.org/wiki/Principal_component_analysis
            %PCA of a multivariate Gaussian distribution centered at (1,3) with a standard deviation of 
            %3 in roughly the (0.878, 0.478) direction and of 1 in the orthogonal direction. 
            %The vectors shown are the eigenvectors of the covariance matrix scaled by the square root of the corresponding eigenvalue, 
            %and shifted so their tails are at the mean.
            %so the matrix of the eigenvectors can rotate the x,y in the new
            %coordinate (i.e. 1,0) to the old coordiate (0.878, 0.478)
            %     X*V=T, V is matrix of the eigenvectors or weights that transform
            %     x,y to m,n in the new coordinate of T and m,n cov is diagonal
            %     matrix of the eighenvalues
            %     V*T'=X' (Or T*V'=X), in the way, the m,n in T can be transform back to x,y in
            %     the original coordinate system. Note that because  X*V=T, T*V'=X

           % VV = V*sqrt(D);               %# scale eigenvectors, VV is also the transformation matrix that transform the 1, 0 vector and 0, 1 vector to the eignevectors
            %VV = sqrt(D); 
            
            %e1=V*sqrt(D)*e+Mu', should we get eS=sqrt(D)*e first then V*eS
            
            %e=sqrt(D)*e
            %e=V*e
            VV = V*sqrt(D);
            e=VV*e;
            
            e = bsxfun(@plus, e, Mu'); %#' project circle back to orig space

            %# plot cov and major/minor axes
            plot(Mu(1),Mu(2),mark,'Color',color);
            hold on;
            plot(e(1,:), e(2,:), 'Color',color);
            %#quiver(Mu(1),Mu(2), VV(1,1),VV(2,1), 'Color','k')
            %#quiver(Mu(1),Mu(2), VV(1,2),VV(2,2), 'Color','k')
       % end
       output.V=V;
       output.D=D;
       output.Mu=Mu;
 
end

