这是一个快速演练。首先,我们创建一个隐藏变量(或“因素”)的矩阵。它有 100 个观测值,并且有两个独立的因素。
>> factors = randn(100, 2);
现在创建一个载荷矩阵。这会将隐藏变量映射到观察到的变量上。假设你观察到的变量有四个特征。那么你的负载矩阵需要是4 x 2
>> loadings = [
1 0
0 1
1 1
1 -1 ];
这告诉您第一个观察到的变量负载在第一个因素上,第二个负载在第二个因素上,第三个变量负载在因素的总和上,第四个变量负载在因素的差异上。
现在创建您的观察结果:
>> observations = factors * loadings' + 0.1 * randn(100,4);
我添加了少量随机噪声来模拟实验误差。pca
现在我们使用统计工具箱中的函数执行 PCA :
>> [coeff, score, latent, tsquared, explained, mu] = pca(observations);
变量score
是主成分分数的数组。这些将通过构造正交,您可以检查 -
>> corr(score)
ans =
1.0000 0.0000 0.0000 0.0000
0.0000 1.0000 0.0000 0.0000
0.0000 0.0000 1.0000 0.0000
0.0000 0.0000 0.0000 1.0000
该组合score * coeff'
将重现您观察的中心版本。mu
在执行 PCA 之前减去平均值。要重现您的原始观察结果,您需要将其重新添加,
>> reconstructed = score * coeff' + repmat(mu, 100, 1);
>> sum((observations - reconstructed).^2)
ans =
1.0e-27 *
0.0311 0.0104 0.0440 0.3378
要获得原始数据的近似值,您可以开始从计算的主成分中删除列。为了了解要删除哪些列,我们检查explained
变量
>> explained
explained =
58.0639
41.6302
0.1693
0.1366
这些条目告诉您每个主成分解释了多少百分比的方差。我们可以清楚地看到,前两个分量比后两个分量更显着(它们解释了它们之间超过 99% 的差异)。使用前两个分量来重建观察结果给出了 rank-2 近似值,
>> approximationRank2 = score(:,1:2) * coeff(:,1:2)' + repmat(mu, 100, 1);
我们现在可以尝试绘制:
>> for k = 1:4
subplot(2, 2, k);
hold on;
grid on
plot(approximationRank2(:, k), observations(:, k), 'x');
plot([-4 4], [-4 4]);
xlim([-4 4]);
ylim([-4 4]);
title(sprintf('Variable %d', k));
end
我们几乎完美地再现了原始观察结果。如果我们想要更粗略的近似,我们可以只使用第一个主成分:
>> approximationRank1 = score(:,1) * coeff(:,1)' + repmat(mu, 100, 1);
并绘制它,
>> for k = 1:4
subplot(2, 2, k);
hold on;
grid on
plot(approximationRank1(:, k), observations(:, k), 'x');
plot([-4 4], [-4 4]);
xlim([-4 4]);
ylim([-4 4]);
title(sprintf('Variable %d', k));
end
这一次的重建不是那么好。那是因为我们故意将我们的数据构建为具有两个因素,而我们只是从其中一个因素中重建它。
请注意,尽管我们构建原始数据的方式与其复制之间存在暗示性的相似性,
>> observations = factors * loadings' + 0.1 * randn(100,4);
>> reconstructed = score * coeff' + repmat(mu, 100, 1);
factors
与之间score
或loadings
与之间不一定有任何对应关系coeff
。PCA 算法对您的数据的构造方式一无所知——它只是试图尽可能多地解释每个连续分量的总方差。
用户@Mari 在评论中询问她如何将重建误差绘制为主成分数量的函数。使用上面的变量explained
非常容易。我将生成一些具有更有趣的因子结构的数据来说明效果 -
>> factors = randn(100, 20);
>> loadings = chol(corr(factors * triu(ones(20))))';
>> observations = factors * loadings' + 0.1 * randn(100, 20);
现在所有的观察都加载在一个重要的共同因素上,而其他因素的重要性逐渐降低。我们可以像以前一样得到 PCA 分解
>> [coeff, score, latent, tsquared, explained, mu] = pca(observations);
并绘制解释方差的百分比如下,
>> cumexplained = cumsum(explained);
cumunexplained = 100 - cumexplained;
plot(1:20, cumunexplained, 'x-');
grid on;
xlabel('Number of factors');
ylabel('Unexplained variance')