Julia中的探索性主成分分析

3
我将尝试用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,如何轻松检索个体坐标?如何轻松显示相关圆等内容?
据我所知,从准机器学习的角度来看,Julia实现了PCA作为“模型”,基本上没有办法轻松地显示所有通常需要在探索性PCA中使用的图形、统计和见解(我的意思是cos²、贡献等)。
是否有其他Julia包更加面向探索性多元分析,能够提供类似于著名的R软件包FactoMineR或ade4的粗略等效物?
谢谢!

1
这里有一个相关的讨论:https://discourse.julialang.org/t/pca-output/22687 - call me Steve
谢谢!这个帖子很有帮助。但是似乎确认了今天在Julia生态系统中与PCA相关的东西确实处于非常基础的状态。 - Philopolis
1个回答

13

首先,您正在使用错误的方向输入数据。正如您可以从文档中看到的那样:

fit(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类型的对象。事实证明,在文档中详细描述了此类型的可用方法:

属性

设M是PCA的实例,d为观察值的维度,p为输出维度(即主子空间的维度)

indim(M)
获取输入维数d,即观测空间的维数。

outdim(M)
获取输出维数p,即主子空间的维数。

mean(M)
获取均值向量(长度为d)。

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

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

principalvars(M)
主成分的方差。

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

tresidualvar(M)
总剩余方差。

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

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

但是,如果您不知道这一点或者在找文档时遇到困难,您可以通过使用methodswith函数获取所有可用方法的完整列表(对于Julia中的任何类型)-- 在本例中特别是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,可以使用它将观测值转换为主成分,如下所示:

= ᵀ(-)

或使用它来从主成分(近似)重构观测数据,如下所示:

̃ = +

这里, 是投影矩阵。

该软件包提供了相应的方法:

transform(M, x)
将观测数据x转换为主成分。

其中,x可以是长度为d的向量或每列代表一个观测的矩阵。

reconstruct(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

鸢尾花数据集的PCA,绘制PC1 vs PC2

最后,如果您想为不同变量添加箭头,那么我们已经在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适用于PCA类型的配方,但是有一个适用于MDS的方法,因此您可以通过简单地编写以下代码获得与上述相同的图:

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

(在这种情况下,当距离度量为原始(即高维)向量空间中的欧几里得距离时,PCA和MDS是等价的)


PCA和MDS在这种情况下是等价的,因为距离度量采用原始向量空间中的欧几里得距离。

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接