PCA 双图 - Matlab

Question

基于答案 Unable to create an array from a table, Coloring The Dots in biPlot Chart我写了下面的代码。 如何更改图例以显示clusters = Tb.class（鸢尾属物种）以及如何对每个组进行凸包？

代码：

clc;
clear;
close all;

Tb = webread('https://datahub.io/machine-learning/iris/r/iris.csv');
clusters = Tb.class;

X = [Tb.sepallength Tb.sepalwidth   Tb.petallength  Tb.petalwidth ];

Z = zscore(X); % Standardized data
[coefs,score] = pca(Z);

vbls = {'sepallength','sepalwidth','petallength','petalwidth'}; 
h=biplot(coefs(:,1:2),'Scores',score(:,1:2),'VarLabels',vbls);

hID = get(h, 'tag'); 
% Isolate handles to scatter points
hPt = h(strcmp(hID,'obsmarker')); 
% Identify cluster groups
grp = findgroups(clusters);    %r2015b or later - leave comment if you need an alternative
grp(isnan(grp)) = max(grp(~isnan(grp)))+1; 
grpID = 1:max(grp); 
% assign colors and legend display name
clrMap = lines(length(unique(grp)));   % using 'lines' colormap
for i = 1:max(grp)
    set(hPt(grp==i), 'Color', clrMap(i,:), 'DisplayName', sprintf('Cluster %d', grpID(i)))
end
% add legend to identify cluster
[~, unqIdx] = unique(grp);
legend(hPt(unqIdx)) 

Answer 1

1) 对于legend，您可以简单地使用集群名称，如下所示：

unqClusters = unique(clusters);
for i = 1:max(grp)
    % Original legend showing "Cluster n"
    %set(hPt(grp==i), 'Color', clrMap(i,:), 'DisplayName', sprintf('Cluster %d', grpID(i)))
    % New legend showing the flower class name
    set(hPt(grp==i), 'Color', clrMap(i,:), 'DisplayName', sprintf('%s', unqClusters{i}))
end

2) 对于凸包部分，可以使用二维正态分布的contours。

这是适用于前两个原始 PC 分数的代码，您需要适应 biplot() 函数执行的标准化。

代码基本上遍历每个集群，计算前两个 PC 值的均值和协方差，并将这些值用作二维正态分布的参数，其 PDF 是在网格网格上计算的。最后选择二维 PDF 的一个级别值进行绘图，这会导致椭圆覆盖集群中的大部分点。

%% Convex hulls around each cluster (ellipses based on 2D normal distribution)
% Execution parameters:
% Number of points used for the meshgrid to generate the Normal PDF
resolution = 100;
% Factor of PC range by which to extend the min and max PC values
% on which the normal PDF is computed
% (adjust this value so that the ellipse around the points is drawn completely)
factor = 2;

% Variables to store the location and scale parameters of each cluster
mu = cell(1,3);
Sigma = cell(1,3);

% Do the calculations for each cluster and plot points and ellipses on a single plot
figure
hold on
for g = grpID
    % Score (PC) values to analyze
    scoregrp = score(grp==g,1:2);
    mu{g} = mean(scoregrp);
    Sigma{g} = cov(scoregrp);

    % PC1 and PC2 values on which to compute the Normal PDF
    pcvalues = cell(1,2);
    for c = 1:length(pcvalues)
        pcrange = range(scoregrp(:,c));
        pcmin = min(scoregrp(:,c)) - factor*pcrange;
        pcmax = max(scoregrp(:,c)) + factor*pcrange;
        pcstep = pcrange / resolution;
        pcvalues{c} = pcmin:pcstep:pcmax;
    end

    % Meshgrid generated from the above PC values
    [XPoints, YPoints] = meshgrid(pcvalues{1}, pcvalues{2});
    XYPoints = [XPoints(:) YPoints(:)];
    npoints = size(XYPoints, 1);

    % Variables used repeatedly in the computation of the normal PDF
    detSigma = det(Sigma{g});
    invSigma = inv(Sigma{g});

    % Normal PDF on the meshgrid
    % Note: the PDF is computed via a loop instead of matrix operation
    % due to memory constraints
    pdfnormal = nan(1, npoints);
    den = sqrt(2*pi*detSigma);
    for i = 1:npoints
        point = XYPoints(i,:) - mu{g};
        pdfnormal(i) = exp( -0.5 * point * invSigma * point' ) / den;
    end

    % Reshape the PDF to 2D matrix so that we can plot the contours
    pdfnormal2D = reshape(pdfnormal, length(pcvalues{1}), length(pcvalues{2}));

    % Value to plot as ellipse where the normal falls at about the value
    % of the standard normal at two standard deviations
    pdfvalue2Std = pdf('norm', 2);

    % Plot!
    plot(scoregrp(:,1), scoregrp(:,2), '.', 'Color', clrMap(g,:))
    contour(XPoints, YPoints, pdfnormal2D, pdfvalue2Std, 'LineWidth', 2, ...
            'Color', clrMap(g,:))
end
xlabel('PC1')
ylabel('PC2')

生成以下图表：