Confidence Ellipsoid from Eigenvalues

Question

Hugo Leevy 2016-8-15

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/299648-confidence-ellipsoid-from-eigenvalues

编辑： Varsha Radhakrishnan 2021-7-14

Hello I am new in the Matlab world.

I am trying to create a 95% Confidence Ellipsoid for a set of data points. For an ellipse I realized it by solving the eigenvalue problem for the covariance matrix of the 2D data points. I derived the information from a script available online:

    % Calculate the eigenvectors and eigenvalues
    covariance = cov(data);
    [eigenvec, eigenval ] = eig(covariance);
    % Get the index of the largest eigenvector
    [largest_eigenvec_ind_c, r] = find(eigenval == max(max(eigenval)));
    largest_eigenvec = eigenvec(:, largest_eigenvec_ind_c);
    % Get the largest eigenvalue
    largest_eigenval = max(max(eigenval));
    % Get the smallest eigenvector and eigenvalue
    if(largest_eigenvec_ind_c == 1)
        smallest_eigenval = max(eigenval(:,2))
        smallest_eigenvec = eigenvec(:,2);
    else
        smallest_eigenval = max(eigenval(:,1))
        smallest_eigenvec = eigenvec(1,:);
    end
    % Calculate the angle between the x-axis and the largest eigenvector
    angle = atan2(largest_eigenvec(2), largest_eigenvec(1));
    % This angle is between -pi and pi.
    % Let's shift it such that the angle is between 0 and 2pi
    if(angle < 0)
        angle = angle + 2*pi;
    end
    % Get the coordinates of the data mean
    avg = mean(data);
    chisquare_val = 2.4477;
    theta_grid = linspace(0,2*pi);
    phi = angle;
    X0=avg(1);
    Y0=avg(2);
    a=chisquare_val*sqrt(largest_eigenval);
    b=chisquare_val*sqrt(smallest_eigenval);

Therefore I could extract the needed parameters to plot the ellipse. Now I tried to convert the same for a ellipsoid. My approach so far:

    covari=cov(data);
    [eigenvec, eigenval ] = eig(covari);
    values=eigenval(eigenval~=0);
    E1=find(values == max(values));
    E3=find(values == min(values));
      if E1 == 1 && E3 == 3
          E2=2;
      end
      if E1 == 1 && E3 == 2
          E2=3;
      end
      if E1 == 2 && E3 == 3
          E2=1;
      end
       if E1 == 2 && E3 == 1
          E2=3;
       end
       if E1 == 3 && E3 == 1
          E2=2;
       end
         if E1 == 3 && E3 == 2
          E2=1;
      end
  first=eigenvec(:,E1);
  second=eigenvec(:,E2);
  third=eigenvec(:,E3);
avg = mean(data);
chisquare_val = 2.4477;
X0=avg(1);
Y0=avg(2);
Z0=avg(3);
a=chisquare_val*sqrt(values(E1));
b=chisquare_val*sqrt(values(E2));
c=chisquare_val*sqrt(values(E3));
[mX,mY,mZ] = ellipsoid(X0,Y0,Z0,a,b,c,30);

Well this creates an ellipsoid with two main problems: It does not cover at all the data points, as happened for instance for the example 3D data points:

data=  
    -1    -8     1
     2   120     2
    -3    11     3
     4     6     4
     5    90     5

And the second problem is the orientation of the Matrix. From PCA I know that the Scores result from a linear combination of the data points with the eigenvector coefficients, thus a projection into the new space defined by the eigenvectors (principal components). Thus I tried the same here and transformed the data points by linear combination.

[r,c]=size(mX);
old1=[mX(:,1) mY(:,1) mZ(:,1)];
newXYZ1=old1*transMat;
newX=[newXYZ1(:,1)];
newY=[newXYZ1(:,2)];
newZ=[newXYZ1(:,3)];
for i=2:c 
oldX=mX(:,i);
oldY=mY(:,i);
oldZ=mZ(:,i);
old=[oldX oldY oldZ];
newXYZ=old*transMat;
newX=[newX newXYZ(:,1)];
newY=[newY newXYZ(:,2)];
newZ=[newZ newXYZ(:,3)];
end

Most likely I did some mistakes, because the obvious overlay of the resulting ellipsoids and the data points are not even touching one another. I am aware that my code is at the moment not really elegant (especially the if statements, but I can optimize it later on, first it should work. I am new here, please be indulgent.

I hope you guys can help me out here. Best, Hugo