Julia 中的探索性 PCA
Exploratory PCA in Julia
我试图了解如何使用包 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 包，它可以提供与 FactoMineR 或 ade4 等著名 R 包的粗略等效？
谢谢！
好的，所以首先你将数据放入错误的方向。从 docs
可以看出
fit(PCA, X; ...)
Perform PCA over the data given in a matrix X. Each column of X is an observation.
您需要将观察作为列，将变量作为行。考虑到我们通常将变量视为列，这可能看起来令人困惑，但在底层线性代数的上下文中它更有意义。因此，为了做到这一点，让我们从以下开始：
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)
如您所见，此 returns 类型为 PCA 的对象。事实证明，这种类型的可用方法在 documentation:
中有很好的描述
Properties
Let M be an instance of PCA, d be the dimension of observations, and p be the output dimension (i.e the dimension of the principal subspace)
indim(M)
Get the input dimension d, i.e the dimension of the observation space.
outdim(M)
Get the output dimension p, i.e the dimension of the principal subspace.
mean(M)
Get the mean vector (of length d).
projection(M)
Get the projection matrix (of size (d, p)). Each column of the projection matrix corresponds to a principal component.
The principal components are arranged in descending order of the corresponding variances.
principalvars(M)
The variances of principal components.
tprincipalvar(M)
The total variance of principal components, which is equal to sum(principalvars(M)).
tresidualvar(M)
The total residual variance.
tvar(M)
The total observation variance, which is equal to tprincipalvar(M) + tresidualvar(M).
principalratio(M)
The ratio of variance preserved in the principal subspace, which is equal to tprincipalvar(M) / tvar(M).
但是如果您不知道或者找不到文档，您可以通过以下方式获得所有可用方法的完整列表（对于 any 类型的 Julia）使用 methodswith 函数——在本例中特别是 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
正如文档告诉我们的那样
each column of the projection matrix corresponds to a principal component
所以 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
展望相同文档的“转换”部分
Transformation and Construction
Given a PCA model M, one can use it to transform observations into principal components, as
= ᵀ ( - )
or use it to reconstruct (approximately) the observations from principal components, as
̃ = +
Here, is the projection matrix.
The package provides methods to do so:
transform(M, x)
Transform observations x into principal components.
Here, x can be either a vector of length d or a matrix where each column is an observation.
reconstruct(M, y)
Approximately reconstruct observations from the principal components given in y.
Here, y can be either a vector of length p or a matrix where each column gives the principal components for an observation.
我们可以看到将此转换应用于我们的数据的方法是
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 在距离度量为原始（即高维）向量中的欧氏距离的情况下等效 space）
我试图了解如何使用包 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 包，它可以提供与 FactoMineR 或 ade4 等著名 R 包的粗略等效？
谢谢！
好的，所以首先你将数据放入错误的方向。从 docs
可以看出fit(PCA, X; ...)
Perform PCA over the data given in a matrix X. Each column of X is an observation.
您需要将观察作为列，将变量作为行。考虑到我们通常将变量视为列，这可能看起来令人困惑，但在底层线性代数的上下文中它更有意义。因此，为了做到这一点，让我们从以下开始：
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)
如您所见，此 returns 类型为 PCA 的对象。事实证明，这种类型的可用方法在 documentation:
Properties
Let M be an instance of PCA, d be the dimension of observations, and p be the output dimension (i.e the dimension of the principal subspace)
indim(M)
Get the input dimension d, i.e the dimension of the observation space.outdim(M)
Get the output dimension p, i.e the dimension of the principal subspace.mean(M)
Get the mean vector (of length d).projection(M)
Get the projection matrix (of size (d, p)). Each column of the projection matrix corresponds to a principal component.The principal components are arranged in descending order of the corresponding variances.
principalvars(M)
The variances of principal components.tprincipalvar(M)
The total variance of principal components, which is equal to sum(principalvars(M)).tresidualvar(M)
The total residual variance.tvar(M)
The total observation variance, which is equal to tprincipalvar(M) + tresidualvar(M).principalratio(M)
The ratio of variance preserved in the principal subspace, which is equal to tprincipalvar(M) / tvar(M).
但是如果您不知道或者找不到文档，您可以通过以下方式获得所有可用方法的完整列表（对于 any 类型的 Julia）使用 methodswith 函数——在本例中特别是 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
正如文档告诉我们的那样
each column of the projection matrix corresponds to a principal component
所以 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
展望相同文档的“转换”部分
Transformation and Construction
Given a PCA model M, one can use it to transform observations into principal components, as
= ᵀ ( - )
or use it to reconstruct (approximately) the observations from principal components, as
̃ = +
Here, is the projection matrix.
The package provides methods to do so:
transform(M, x)
Transform observations x into principal components.Here, x can be either a vector of length d or a matrix where each column is an observation.
reconstruct(M, y)
Approximately reconstruct observations from the principal components given in y.Here, y can be either a vector of length p or a matrix where each column gives the principal components for an observation.
我们可以看到将此转换应用于我们的数据的方法是
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 在距离度量为原始（即高维）向量中的欧氏距离的情况下等效 space）