这是一个快速的演练.首先我们创建一个隐藏变量的矩阵.它有100个观察,有两个特点.
>> Y = randn(100,2);
现在创建一个加载矩阵.这将把隐藏的变量映射到观察到的变量上.说你观察到的变量有四个特征.那么你的加载矩阵需要是2×4
>> W = [1 1; 1 -1; 2 1; 2 -1]';
这告诉你,观察的第一部分是两个因素的总和,第二个是两个因素的差异,第三个和第四个是各种其他组合.
现在创建你的观察:
>> X = Y * W + 0.1 * randn(100,4);
我添加了少量随机噪声来模拟实验误差.现在我们使用stats工具箱中的princomp函数执行PCA:
>> [w pc ev] = princomp(X);
你应该如何解释这些?那么,pc是主要组件的矩阵.它应该拉出与原始Y变量非常接近的因子.你可以查看这个:
>> corr(pc(:,1),Y(:,1)); % returns -0.9981
>> corr(pc(:,2),Y(:,2)); % returns 0.9830
组合pc * w’将重新创建您的原始数据,减去其平均值.在执行PCA之前总是减去平均值.所以要获取我们做的原始数据
>> mu = mean(X);
>> xhat = bsxfun(@minus,X,mu); % subtract the mean
>> norm(pc * w' - xhat);
ans =
2.3385e-014
因为w是正交的,所以你也有Xhat * w = pc或者示意性地(即这个代码不会执行)
(X - mu) * w = pc <=> X = mu + pc * w'
要获得原始数据的近似值,您可以从计算的主要组件开始删除列.要了解哪些列要放弃,我们检查ev变量
>> ev
ev =
11.1323
3.0812
0.0116
0.0068
我们可以清楚地看到,前两个因素比第二个因素更重要.所以我们来试试看
>> Xapprox = pc(:,1:2) * w(:,1:2)';
>> Xapprox = bsxfun(@plus,mu,Xapprox); % add the mean back in
我们现在可以尝试绘图:
>> plot(Xapprox(:,1),X(:,1),'.'); hold on; plot([-4 4],[-4 4])
>> xlabel('Approximation'); ylabel('Actual value'); grid on;
它看起来像一个非常合理的近似值.它高估了一点,但我们不能都是完美的.
如果我们想要一个更粗略的近似,我们可以使用第一个主成分:
>> Xapprox = pc(:,1) * w(:,1)';
>> Xapprox = bsxfun(@plus,mu,Xapprox);
>> plot(Xapprox(:,1),X(:,1),'.'); hold on; plot([-4 4],[-4 4])
>> xlabel('Approximation'); ylabel('Actual value'); grid on;
这次重建不是很好.那是因为我们故意构造我们的数据有两个因素,我们只是从其中一个重建它们.
最后,您可能希望看到每个因素解释了多少差异.您可以使用ev变量来执行此操作:
>> 100*ev/sum(ev)
ans =
78.2204
21.6499
0.0818
0.0479
因此,第一个组件解释了78%的方差,下一个组件解释了大约22%,最后两个组件解释了其余部分.