修复MATLAB中带95%置信椭球的3D PCA绘图代码错误

问题背景 #

需要绘制包含前三个主成分的3D PCA散点图，并为每个类别添加95%置信椭球，但MATLAB运行代码时抛出reshape错误：

Error using reshape Number of elements must not change. Use [] as one of the size inputs to automatically calculate the appropriate size for that dimension.
Error in test_3Dpcav2>drawEllipsoid (line 83)
points = reshape(points, size(x,1), size(x,2), size(x,3));
Error in test_3Dpcav2 (line 60) drawEllipsoid(mean1, cov1, confLevel, 'r')

错误原因 #

1. Reshape维度不匹配：sphere()默认返回21×21的二维矩阵，size(x,3)结果为1，但处理后的points是441×3的矩阵（总元素1323个），而reshape(points,21,21,1)仅需441个元素，元素总数不匹配。
2. 置信区间计算错误：误用F分布逆函数finv，多元正态分布的置信椭球临界值应使用卡方分布逆函数chi2inv，自由度等于特征维度（此处为3）。
3. 变量名冲突：使用MATLAB内置函数名mean作为变量，可能引发潜在逻辑冲突。

修复后的完整代码 #

% Input data
InputData = [7.72 6.73 3.33 0.12 0.06 -0.31 1;
             8.92 8.22 4.56 0.06 0.72 0.01 1;
             2.11 1.59 0.7 0.12 0.67 -0.02 1;
             79.42 92.88 75.03 64.98 67.8 55.93 2;
             100.82 112.5 94.7 73.16 74.22 64.15 2;
             69.43 83.4 67.46 80.78 81.36 69.38 2;
             25.25 23.86 17.95 0.12 0.44 0.14 3;
             39.97 36.9 30.06 0.13 0.82 0.35 3;
             31.99 29.48 23.68 -0.12 0.53 -0.1 3;
             129.93 126.54 113.51 1.09 0.17 0.38 4;
             129.91 127.4 114.6 1.32 1.05 1.19 4;
             110.38 108.7 96.3 0.41 1.25 0.87 4];

% Extract the features (first 6 columns) and labels (last column) from the data
X = InputData(:,1:6);
labels = InputData(:,7);

% Perform PCA on the data
[coeff,score,latent,~,explained] = pca(X);

% Get the first three principal components
PC1 = score(:,1);
PC2 = score(:,2);
PC3 = score(:,3);

% Calculate the means for each label
mean1 = mean(score(labels==1,1:3));
mean2 = mean(score(labels==2,1:3));
mean3 = mean(score(labels==3,1:3));
mean4 = mean(score(labels==4,1:3));

% Calculate the covariance matrices for each label
cov1 = cov(score(labels==1,1:3));
cov2 = cov(score(labels==2,1:3));
cov3 = cov(score(labels==3,1:3));
cov4 = cov(score(labels==4,1:3));

% Set the confidence level for the ellipsoids
confLevel = 0.95;

% Create a 3D scatter plot with different colors for each label
figure
scatter3(PC1(labels==1), PC2(labels==1), PC3(labels==1), 'r', 'filled')
hold on
scatter3(PC1(labels==2), PC2(labels==2), PC3(labels==2), 'g', 'filled')
scatter3(PC1(labels==3), PC2(labels==3), PC3(labels==3), 'b', 'filled')
scatter3(PC1(labels==4), PC2(labels==4), PC3(labels==4), 'm', 'filled')
title('3D Scatter Plot with 95% Confidence Ellipsoids')
xlabel('PC1')
ylabel('PC2')
zlabel('PC3')
legend('Label 1', 'Label 2', 'Label 3', 'Label 4')

% Create the ellipsoids for each label
drawEllipsoid(mean1, cov1, confLevel, 'r')
drawEllipsoid(mean2, cov2, confLevel, 'g')
drawEllipsoid(mean3, cov3, confLevel, 'b')
drawEllipsoid(mean4, cov4, confLevel, 'm')

% Function to draw ellipsoids
function drawEllipsoid(mu, covMatrix, confLevel, color)
    % Get the eigenvalues and eigenvectors of the covariance matrix
    [V,D] = eig(covMatrix);

    % Calculate the radii of the ellipsoid using chi-squared distribution
    chi2_val = chi2inv(confLevel, size(covMatrix,1));
    radii = sqrt(diag(D) * chi2_val);

    % Generate points on a unit sphere
    [x,y,z] = sphere;

    % Transform the points using the eigenvectors and radii
    points = [x(:) y(:) z(:)] * diag(radii) * V';

    % Add the mean to the points
    points = points + mu;

    % Reshape the points into 3D matrix matching sphere's dimensions
    points = reshape(points, size(x,1), size(x,2), 3);

    % Plot the ellipsoid
    h = surf(points(:,:,1), points(:,:,2), points(:,:,3));
    set(h, 'FaceColor', color, 'FaceAlpha', 0.1, 'EdgeAlpha', 0.1)
end

关键修复点 #

• 将reshape的第三个维度改为3，匹配points的3列特征，保证元素总数一致（21×21×3=1323，与441×3元素数相同）。
• 替换finv为chi2inv，使用卡方分布计算置信椭球的临界值，符合多元正态分布的统计逻辑。
• 将函数内的变量名mean改为mu，避免覆盖MATLAB内置函数。

内容的提问来源于stack exchange，提问作者STEMQs