Julia 中的探索性 PCA答案

【问题标题】：Exploratory PCA in JuliaJulia 中的探索性 PCA
【发布时间】：2021-06-20 08:03:24
【问题描述】：

我尝试了解如何使用 MultivariateStats.jl 包在 Julia 中执行简单的探索性 PCA。

例如，在 R 中，可以执行以下操作：

library(FactoMineR)
data(iris)

## PCA model:
res_pca <- PCA(iris, quali.sup = 5, graph = FALSE)
## Retrieve coordinates of individuals on all 4 PCs:
res_pca$ind$coord
## Simple plots:
plot(res_pca, choix = "ind", habillage = 5)
plot(res_pca, choix = "var")

我无法在 Julia 中获得任何等效的这些非常基本的操作（我的意思是使用 MultivariateStats.jl 中的“本机”函数）。让我们开始吧：

using MultivariateStats
using RDatasets

iris = dataset("datasets", "iris")

## PCA model:
M = fit(PCA, Array(iris[:, 1:4]); pratio=1, maxoutdim=4)

一旦获得 PCA“模型”M，我如何轻松检索个人坐标？如何轻松显示相关圈等？

据我了解，从准机器学习的角度来看，PCA 在 Julia 中实现为“模型”，并且基本上无法轻松显示探索性 PCA 中通常需要的所有常用图表、统计数据和见解？（我的意思是 cos²、贡献等）

是否有任何其他 Julia 包更倾向于探索性多元分析，可以提供与著名 R 包（如 FactoMineR 或 ade4）大致相当的功能？

谢谢！

【问题讨论】：

这里有相关讨论：discourse.julialang.org/t/pca-output/22687
谢谢！这个线程很有帮助。但这似乎证实了 Julia 生态系统中与 PCA 相关的东西确实处于非常初级的状态。

标签： julia pca

【解决方案1】：

好的，首先您将数据放入错误的方向。从docs可以看到

拟合（PCA，X；...）
对矩阵 X 中给定的数据执行 PCA。X 的每一列都是一个观察值。

您需要将观察作为列，将变量作为行。鉴于我们通常将变量视为列，这可能看起来令人困惑，但在底层线性代数的上下文中它更有意义。因此，为了做到这一点，让我们从以下开始：

using MultivariateStats, Statistics
using RDatasets: dataset

iris = dataset("datasets", "iris")
iris_matrix = Array(iris[:, 1:4])'

## PCA model:
M = fit(PCA, iris_matrix; pratio=1, maxoutdim=4)

如您所见，这会返回一个PCA 类型的对象。事实证明，这种类型的可用方法在documentation 中有很好的描述：

属性

令 M 为 PCA 的一个实例，d 为观测的维度，p 为输出维度（即主子空间的维度）

indim(M)
获取输入维度d，即观察空间的维度。

暗淡(M)
获取输出维度p，即主子空间的维度。

平均（米）
获取平均向量（长度为 d）。

投影（米）
获取投影矩阵（大小为 (d, p)）。投影矩阵的每一列对应一个主成分。

主成分按对应方差的降序排列。

校长(M)
主成分的方差。

tprincipalvar(M)
主成分的总方差，等于sum(principalvars(M))。

tresidualvar(M)
总残差。

tvar(M)
总观测方差，等于 tprincipalvar(M) + tresidualvar(M)。

本金比(M)
主子空间中保留的方差比，等于 tprincipalvar(M) / tvar(M)。

但如果您不知道或找不到文档，您可以使用 methodswith 获取所有可用方法的完整列表（如 Julia 中的 any 类型）函数——在这种情况下特别是methodswith(PCA)。

其中的关键是projection，它为我们提供了实际的投影矩阵：

julia> proj = projection(M)
4×4 Matrix{Float64}:
 -0.361387    0.656589  -0.58203     0.315487
  0.0845225   0.730161   0.597911   -0.319723
 -0.856671   -0.173373   0.0762361  -0.479839
 -0.358289   -0.075481   0.545831    0.753657

正如文档告诉我们的那样

投影矩阵的每一列对应一个主成分

所以proj[:,1] 包含 PC1 的权重/“负载”，proj[:,2] 包含 PC2 的负载等。

由于您询问了贡献，我猜您的意思是每个主成分对解释总方差的贡献，文档告诉我们您可以通过 principalvars 得到它：

julia> principalvars(M)
4-element Vector{Float64}:
 4.2282417060348605
 0.24267074792863352
 0.07820950004291898
 0.023835092973449976

所以 PC1 确实在这里完成了大部分繁重的工作。或者，如果您更喜欢每个组件解释的总方差百分比，那么：

julia> principalvars(M) ./ tvar(M) * 100
4-element Vector{Float64}:
 92.46187232017269
  5.306648311706788
  1.7102609807929683
  0.5212183873275495

展望同一文档的“转换”部分

改造与建设

给定一个 PCA 模型 M，可以使用它将观察值转换为主成分，如

? = ?ᵀ (? - ?)

或使用它来重建（近似）主成分的观察结果，如

?̃ = ? ? + ?

这里，?是投影矩阵。

该包提供了这样做的方法：

变换(M, x)
将观测值 x 转换为主成分。

这里，x 可以是一个长度为 d 的向量，也可以是一个矩阵，其中每列都是一个观察值。

重构(M, y)
根据 y 中给出的主成分近似重建观测值。

在这里，y 可以是长度为 p 的向量，也可以是矩阵，其中每列给出观察的主成分。

我们可以看到将这种转换应用于我们的数据的方式是

iris_transformed = transform(M, iris_matrix)

或者如果你想自己做线性代数

iris_transformed = projection(M)' * (iris_matrix .- mean(M))

然后，一旦我们获得了转换后的数据，我们就可以将前两个主成分相互对比

using Plots
h = plot(iris_transformed[1,:], iris_transformed[2,:], seriestype=:scatter, label="")
plot!(xlabel="PC1", ylabel="PC2", framestyle=:box) # A few formatting options

最后，如果你想为不同的变量添加箭头，事实证明我们已经在projection(M) 中拥有了我们需要的一切，我们已经将其存储为proj

for i=1:4; plot!([0,proj[i,1]], [0,proj[i,2]], arrow=true, label=names(iris)[i], legend=:bottomleft); end
display(h)

编辑：我可能还应该提到，虽然（还）没有用于 PCA 类型的 StatsPlots 配方，但有一个用于 MDS 的配方，因此在这种情况下，您可以通过简单编写获得与上述等效的图

M = fit(MDS, iris_matrix; distances=false)
using StatsPlots
plot(M)

（对于此处距离度量是原始（即高维）向量空间中的欧几里得距离的情况，PCA 和 MDS 是等价的）

【讨论】：