# A question ahout PCA

**URL:** <https://discourse.julialang.org/t/a-question-ahout-pca/42100>\
**Category:** Optimization (Mathematical)\
**Tags:** question\
**Created:** [June 26, 2020, 1:21pm UTC](https://discourse.julialang.org/t/a-question-ahout-pca/42100 "2020-06-26T13:21:29Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![abcde](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abcde/32/16180_2.png) [@abcde](https://discourse.julialang.org/u/abcde)\
**Post date:** [June 26, 2020, 1:21pm UTC](https://discourse.julialang.org/t/a-question-ahout-pca/42100/1 "2020-06-26T13:21:29Z")

</div>

After applying PCA algorithm, how to display its principal component?thanks!!

 ![image](https://global.discourse-cdn.com/julialang/original/3X/1/c/1c41bfee2cb4be97b9261cff5d45f8b4371c57ce.png)

---

<div class="post-metadata">

**Author:** ![nilshg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nilshg/32/2283_2.png) [@nilshg](https://discourse.julialang.org/u/nilshg)\
**Post date:** [June 26, 2020, 1:45pm UTC](https://discourse.julialang.org/t/a-question-ahout-pca/42100/2 "2020-06-26T13:45:37Z")

</div>

There have been quite a few discussions on this, use the search function on discourse to find e.g.

> [@PCA Output?](https://discourse.julialang.org/t/pca-output/22687/24):
>
> To support this request: …so better check that we don’t choose an approach/terminology which is only familiar to psychologists. The use of the function from the psych package is widespread and I would say it is common knowledge among R users. …ike factor analysis. Factor analysis is also affected. The output of the model estimate is very similar.

---

<div class="post-metadata">

**Author:** ![eliassno](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/eliassno/32/18917_2.png) [@eliassno](https://discourse.julialang.org/u/eliassno)\
**Post date:** [June 26, 2020, 1:48pm UTC](https://discourse.julialang.org/t/a-question-ahout-pca/42100/3 "2020-06-26T13:48:17Z")

</div>

This was already answered by @nilshg in the post mentioned above.

* * *

You need the projection matrix, which you can get with `projection(M)`. You can read more about the [properties in the docs](https://multivariatestatsjl.readthedocs.io/en/latest/pca.html#properties).

```julia
using MultivariateStats
X = rand(4,5)
M = fit(PCA, X; maxoutdim=2)
proj = projection(M)[:, 1]

```

---

<div class="post-metadata">

**Author:** ![abcde](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abcde/32/16180_2.png) [@abcde](https://discourse.julialang.org/u/abcde)\
**Post date:** [June 26, 2020, 2:15pm UTC](https://discourse.julialang.org/t/a-question-ahout-pca/42100/4 "2020-06-26T14:15:43Z")

</div>

thanks，but it does not tell me which is the extracted principal component?it only several numbers.

---

<div class="post-metadata">

**Author:** ![abcde](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abcde/32/16180_2.png) [@abcde](https://discourse.julialang.org/u/abcde)\
**Post date:** [June 26, 2020, 2:17pm UTC](https://discourse.julialang.org/t/a-question-ahout-pca/42100/5 "2020-06-26T14:17:30Z")

</div>

Because I want to plot these principal components

---

<div class="post-metadata">

**Author:** ![Gunter\_Faes](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gunter_faes/32/12580_2.png) [@Gunter\_Faes](https://discourse.julialang.org/u/Gunter_Faes)\
**Post date:** [June 26, 2020, 3:50pm UTC](https://discourse.julialang.org/t/a-question-ahout-pca/42100/6 "2020-06-26T15:50:32Z")

</div>

@abcde, I understand you, and I got into this a while back. I was not completely satisfied and so I addressed R about Julia. Here is just a code snippet:

```julia
using RCall
...
R"library(psych)"
# Estimate the number of HK's using the R function scree:
HK_R = R"scree($PCA_Daten)"

...

PCA_Modell = R"principal($PCA_Daten, $PCA_HK)"
PCA_scores = R"principal($PCA_Daten, $PCA_HK)$scores"
show(PCA_Modell)
PCA_scores_M = convert(Matrix, PCA_scores)
...

```

And you will get something like this:

```julia
RObject{VecSxp}
Principal Components Analysis
Call: principal(r = `#JL`$PCA_Daten, nfactors = 4)
Standardized loadings (pattern matrix) based upon correla-tion matrix
RC1 RC2 RC3 RC4 h2 u2 com
1 0.91 -0.13 0.07 0.02 0.85 0.15 1.1
2 -0.34 -0.02 -0.72 -0.13 0.64 0.36 1.5
3 0.75 0.04 0.44 0.23 0.81 0.19 1.9
4 0.35 0.46 -0.03 -0.25 0.39 0.61 2.5
5 0.13 0.00 -0.27 0.80 0.73 0.27 1.3
6 -0.14 0.88 0.03 0.07 0.80 0.20 1.1
7 -0.02 0.88 -0.11 0.09 0.79 0.21 1.1
8 0.78 0.10 -0.44 0.01 0.80 0.20 1.6
9 -0.73 0.01 -0.01 -0.25 0.60 0.40 1.2
10 0.15 0.06 0.27 0.76 0.68 0.32 1.4
11 -0.23 -0.11 0.78 -0.13 0.69 0.31 1.3
RC1 RC2 RC3 RC4
SS loadings 2.87 1.80 1.67 1.44
Proportion Var 0.26 0.16 0.15 0.13
Cumulative Var 0.26 0.42 0.58 0.71
Proportion Explained 0.37 0.23 0.21 0.19
Cumulative Proportion 0.37 0.60 0.81 1.00
Mean item complexity = 1.4
Test of the hypothesis that 4 components are sufficient.
The root mean square of the residuals (RMSR) is 0.09
with the empirical chi square 1473.01 with prob < 3.2e-303
Fit based upon off diagonal values = 0.89

```

But I have to admit that I have not tried out what Julia now delivers as an output. So this is my quick answer…! 😉

---

<div class="post-metadata">

**Author:** ![eliassno](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/eliassno/32/18917_2.png) [@eliassno](https://discourse.julialang.org/u/eliassno)\
**Post date:** [June 26, 2020, 4:41pm UTC](https://discourse.julialang.org/t/a-question-ahout-pca/42100/7 "2020-06-26T16:41:51Z")

</div>

Fair enough. Each column of `proj` gives the corresponding PCA axis. I assumed you only wanted the first axis, since you didn’t provide the data you wanted to transform.

---

<div class="post-metadata">

**Author:** ![eliassno](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/eliassno/32/18917_2.png) [@eliassno](https://discourse.julialang.org/u/eliassno)\
**Post date:** [June 26, 2020, 4:44pm UTC](https://discourse.julialang.org/t/a-question-ahout-pca/42100/8 "2020-06-26T16:44:17Z")

</div>

To visualize your data in 2D with PCA, you can `transform` your data `X`.

```julia
using PyCall, MultivariateStats, Plots

data = pyimport("sklearn.datasets")

# Generate labeled data with 3 classes
X, y = data.make_blobs(n_samples=1000, n_features = 10, centers = 3, cluster_std = 3, random_state=0)

# Columns correspond to observations
X = Array(X')

# Fit PCA model in 2d, for visualization
M = fit(PCA, X; maxoutdim=2)
             
# Transform data
x = transform(M, X)

# Plot dims 1 and 2 w/o PCA
p1 = scatter(X[1,:], X[2,:], c=y, xlabel="Dim. 1", ylabel="Dim. 2", title="Original")

# Plot dims 1 and 2 with PCA
p2 = scatter(x[1,:], x[2,:], xlabel="PC1", ylabel="PC2", c=y, title="PCA")

# Compare plots
plot(p1, p2, layout=2, legend=false)

```

![pca](https://global.discourse-cdn.com/julialang/original/3X/1/9/19b1f4e2ef68c2cc1f0285f98bb9858a517d72f4.png)

---

<div class="post-metadata">

**Author:** ![sylvaticus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sylvaticus/32/203883_2.png) [@sylvaticus](https://discourse.julialang.org/u/sylvaticus)\
**Post date:** [June 26, 2020, 5:19pm UTC](https://discourse.julialang.org/t/a-question-ahout-pca/42100/9 "2020-06-26T17:19:16Z")

</div>

In [BetaML](https://github.com/sylvaticus/BetaML.jl) you have what I believe is one of the simplest approach (disclaimer: I am the author):

“”’

pca(X;K,error)

Perform Principal Component Analysis returning the matrix reprojected among the dimensions of maximum variance.

# Parameters:

- `X` : The (N,D) data to reproject

- `K` : The number of dimensions to maintain (with K\<=D) [def: `nothing`]

- `error`: The maximum approximation error that we are willing to tollerate [def: `0.05`]

# Return:

- A named tuple with:

- `X`: The reprojected (NxK) matrix

- `K`: The dimensions retieved

- `error`: The actual proportion of variance not explained in the reprojected dimensions

- `P`: The (D,K) matrix of the eigenvectors associated to the K-largest eigenvalues used to reproject the data matrix

Note that if `K` is indicated, the parameter `error` has no effect.

“”"

---

<div class="post-metadata">

**Author:** ![colintbowers](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/colintbowers/32/8033_2.png) [@colintbowers](https://discourse.julialang.org/u/colintbowers)\
**Post date:** [June 29, 2020, 1:04am UTC](https://discourse.julialang.org/t/a-question-ahout-pca/42100/10 "2020-06-29T01:04:23Z")

</div>

If the way the package works is bothering you, it’s only a few lines of code to get them yourself using `LinearAlgebra`. For example, in most cases I’m familiar with, users want the components ordered from largest eignvalue to smallest. Given data matrix in `x` with rows as observations and columns as variables/dimensions:

```julia
    x = x .- mean(x, dims=1) # Centre the data
    xcov = Symmetric(cov(x)) #Compute the covariance matrix
    ef = eigen(xcov) #Get eigenvalues and eigenvectors
    eigvalvec = reverse(ef.values, dims=1) #Orient eigenvalues largest to smallest
    eigvecmat = reverse(ef.vectors, dims=2) #Orient eigenvectors the same
    eigvalcs = cumsum(eigvalvec) #Cumulative sum of eigenvalues
    vratio = eigvalcs ./ eigvalcs[end] #Ratios of variance explained
    pc = x * eigvecmat #Principal components ordered largest eigenvalue to smallest

```

You didn’t actually ask for the variance ratios but again they are trivial to compute and often useful so I thought I’d include them.

Cheers,

Colin