高校生向けの主成分分析教材.
参考テキスト 統計科学研究所 様.　Rによる主成分分析.

1 分散共分散行列.

1変量確率変数の分散はイメージしやすい。分散が大きければ広がりが大きい。
２変量確率変数の分散共分散行列から分布の状態をイメージできるだろうか。
固有値・固有ベクトルとの関係も含めてシミュレーションしてみたい。
2変量正規分布に対して主成分分析を適用する。

共分散とはX,Yの一致度である。
共分散行列とは？(定義、意味、分布の推定).

1.1 相関のある2変量

簡単のため、2変量正規乱数を使用する。
算出するのは、

  • 分散共分散行列.
  • 固有値.
  • 固有ベクトル.

散布図と固有ベクトルの関係をつかもう。

library(MASS)
mu <- c(0, 0)
Sigma <- matrix(c(1, 0.7, 0.7, 1), 2, 2)
data <- mvrnorm(50, mu, Sigma)
d <- data.frame(X=data[,1],Y=data[,2])
cov(d)#分散共分散行列
##           X         Y
## X 1.3230098 0.9743596
## Y 0.9743596 1.0634820
eigen_d <- eigen(cov(d))#固有値
v <- eigen_d$vectors
e <- eigen_d$values;e
## [1] 2.1762084 0.2102834
plot(d,asp=1)
arrows(x0=0, y0=0, x1=e[1]*v[1,1],  y1=e[1]*v[2,1], code=2, col="red")
arrows(x0=0, y0=0, x1=v[1,2],  y1=v[2,2], code=2, col="blue")

分布は、
　固有ベクトル (1,1) 方向に 1+r 倍.
　固有ベクトル (1,-1) 方向に 1-r 倍.
相関係数rが大きくなれば楕円の短軸は圧縮される。

1.2 prcomp.

library("psych")
#(describe(d))#基本統計情報
par(family = "HiraKakuProN-W3") #日本語フォントの指定
result <- prcomp(d, scale=T) #scale=Tは相関行列から分析
summary(result)
## Importance of components:
##                           PC1     PC2
## Standard deviation     1.3496 0.42257
## Proportion of Variance 0.9107 0.08928
## Cumulative Proportion  0.9107 1.00000
biplot(result)

biplot は主成分得点と主成分負荷量を同時にプロットしたもので、座標軸に２種類の目盛がある。便利であるがテキストをグラフに記入するため表示にズレがあるので、次のように plotで描くこともできる。

データ科学便覧参照.

  • 主成分スコア $x，標準偏差 $sdev (固有値の正の平方根).
  • $rotation (各主成分軸の固有ベクトル).
  • $center (解析に用いた各項目の値の平均).

次元縮約だけなら，$x を考えれば良い．

#str(result)
#result$x[,1]
#result$x[,2]
result$rotation
##          PC1        PC2
## X -0.7071068 -0.7071068
## Y -0.7071068  0.7071068
rr <- result$rotation#固有ベクトル
plot(d,asp=1,type="n")
text(result$x[,1],result$x[,2],1:50) 
abline(a=0,b=rr[2]/rr[1],col="blue")
abline(a=0,b=rr[4]/rr[3],col="red")

#arrows(0,0,rr[1],rr[2],length=0.05,col="red")

$rotation は2つの固有ベクトル方向を示してくれる。
第１主成分方向（青）　　第2主成分方向（薄青）

主成分負荷量は主成分との相関係数であり、rotationの行をベクトル表示したものである。主成分が元の変数に影響する度合いを図から読み取ることができる。
因子負荷量をプロットすることにより、主成分に寄与している因子を視覚的に捉えることができる。

plot(result$x[,1],result$x[,2],type="n",asp=1) 
text(result$x[,1],result$x[,2],1:50) 
abline(a=0,b=rr[2]/rr[1],col="blue")
abline(a=0,b=rr[4]/rr[3],col="red")

主成分の計算を再現する。
直交行列result$rotation の１列目が第１主成分を求める式の係数。
dを標準化したbの１行目が標準化された元のデータ。
内積を計算すると第１主成分となる。

result$x[1,1:2]#主成分
##        PC1        PC2 
## -1.5933630 -0.2982577
a<-result$rotation;a
##          PC1        PC2
## X -0.7071068 -0.7071068
## Y -0.7071068  0.7071068
b<-scale(d)
a[1:2,1]
##          X          Y 
## -0.7071068 -0.7071068
b[1,1:2]
##         X         Y 
## 1.3375779 0.9157777
a[1:2,1]*b[1,1:2]
##          X          Y 
## -0.9458104 -0.6475526
sum(a[1:2,1]*b[1,1:2])
## [1] -1.593363

1.3 因子負荷量.

MARGIN=2 は行の操作を表す。（MARGIN=1は列操作）

result$sdev
## [1] 1.3496051 0.4225707
fc.l <- sweep(result$rotation, MARGIN=2, result$sdev, FUN="*") 
result$rotation
##          PC1        PC2
## X -0.7071068 -0.7071068
## Y -0.7071068  0.7071068
z<-fc.l#主成分負荷行列
z
##          PC1        PC2
## X -0.9543149 -0.2988026
## Y -0.9543149  0.2988026
z[,1]#第１主成分の主成分負荷量
##          X          Y 
## -0.9543149 -0.9543149
sum(z[,1]^2)#大きさは第１固有値
## [1] 1.821434
par(family = "HiraKakuProN-W3") #日本語フォントの指定
subject <- c("X", "Y") ;subject
## [1] "X" "Y"
plot(fc.l[,1], pch=subject, ylim=c(-1,1), main="PC1")

plot(fc.l[,2], pch=subject, ylim=c(-1,1), main="PC2")

#主成分との相関係数

1.4 バリマックス回転の主成分分析.

library(psych)
dat <- d
fit <- principal(dat, nfactors=2, rotate=TRUE)#バリマックス回転
## Specified rotation not found, rotate='none' used
fit
## Principal Components Analysis
## Call: principal(r = dat, nfactors = 2, rotate = TRUE)
## Standardized loadings (pattern matrix) based upon correlation matrix
##    PC1  PC2 h2       u2 com
## X 0.95 -0.3  1 -1.1e-16 1.2
## Y 0.95  0.3  1  0.0e+00 1.2
## 
##                        PC1  PC2
## SS loadings           1.82 0.18
## Proportion Var        0.91 0.09
## Cumulative Var        0.91 1.00
## Proportion Explained  0.91 0.09
## Cumulative Proportion 0.91 1.00
## 
## Mean item complexity =  1.2
## Test of the hypothesis that 2 components are sufficient.
## 
## The root mean square of the residuals (RMSR) is  0 
##  with the empirical chi square  0  with prob <  NA 
## 
## Fit based upon off diagonal values = 1
summary(fit)
## 
## Factor analysis with Call: principal(r = dat, nfactors = 2, rotate = TRUE)
## 
## Test of the hypothesis that 2 factors are sufficient.
## The degrees of freedom for the model is -2  and the objective function was  0 
## The number of observations was  50  with Chi Square =  0  with prob <  NA 
## 
## The root mean square of the residuals (RMSA) is  0
biplot(fit)

library(psych)
dat <- d
fit <- principal(dat, nfactors=1, rotate=TRUE)#バリマックス回転
fit
## Principal Components Analysis
## Call: principal(r = dat, nfactors = 1, rotate = TRUE)
## Standardized loadings (pattern matrix) based upon correlation matrix
##    PC1   h2    u2 com
## X 0.95 0.91 0.089   1
## Y 0.95 0.91 0.089   1
## 
##                 PC1
## SS loadings    1.82
## Proportion Var 0.91
## 
## Mean item complexity =  1
## Test of the hypothesis that 1 component is sufficient.
## 
## The root mean square of the residuals (RMSR) is  0.09 
##  with the empirical chi square  0.8  with prob <  NA 
## 
## Fit based upon off diagonal values = 0.99
summary(fit)
## 
## Factor analysis with Call: principal(r = dat, nfactors = 1, rotate = TRUE)
## 
## Test of the hypothesis that 1 factor is sufficient.
## The degrees of freedom for the model is -1  and the objective function was  0.31 
## The number of observations was  50  with Chi Square =  14.42  with prob <  NA 
## 
## The root mean square of the residuals (RMSA) is  0.09
biplot(fit)

2 3変量.

#3D表示
options(rgl.printRglwidget = TRUE) 
#WebGLに変換してRStudio, VSCode, ウェブブラウザなど出力
library(rgl)

library(MASS)
mu <- c(0,0,0)
Sigma <- matrix(c(1, 0.7,0.6,0.7,1,0.5,0.6,0.5, 1), 3, 3)
data <- mvrnorm(100, mu, Sigma)
d <- data.frame(X=data[,1],Y=data[,2],Z=data[,3])
cov(d)#分散共分散行列
##           X         Y         Z
## X 0.8261667 0.5483191 0.4519298
## Y 0.5483191 0.9852074 0.5141435
## Z 0.4519298 0.5141435 0.8682949
attach(d)
#plot3d(X, Y, Z, type = "s", col = "blue", size = 1)

eigen_d <- eigen(cov(d))#固有値
v <- eigen_d$vectors;v[1,1]
## [1] 0.5482452
e <- eigen_d$values;e[1]
## [1] 1.909684
plot3d(X, Y, Z, type = "s", col = "gray", size = 1)
arrow3d(p0 = c(0, 0, 0), p1 = c(e[1]*v[1,1],e[1]*v[2,1],e[1]*v[3,1]), col = "red")
arrow3d(p0 = c(0, 0, 0), p1 = c(2*e[2]*v[1,2],2*e[2]*v[2,2],2*e[2]*v[3,2]), col = "blue")
arrow3d(p0 = c(0, 0, 0), p1 = c(3*e[3]*v[1,3],3*e[3]*v[2,3],3*e[3]*v[3,3]), col = "green")
library("psych")
#(describe(d))#基本統計情報
par(family = "HiraKakuProN-W3") #日本語フォントの指定
result <- prcomp(d, scale=T) #scale=Tは相関行列から分析
summary(result)
## Importance of components:
##                           PC1    PC2    PC3
## Standard deviation     1.4602 0.6912 0.6247
## Proportion of Variance 0.7107 0.1592 0.1301
## Cumulative Proportion  0.7107 0.8699 1.0000
biplot(result)

plot(result$x[,1],result$x[,2],type="n",asp=1) 
text(result$x[,1],result$x[,2],1:100) 

library(psych)
dat <- d
fit <- principal(dat, nfactors=2, rotate=TRUE)#バリマックス回転
## Specified rotation not found, rotate='none' used
fit
## Principal Components Analysis
## Call: principal(r = dat, nfactors = 2, rotate = TRUE)
## Standardized loadings (pattern matrix) based upon correlation matrix
##    PC1   PC2   h2     u2 com
## X 0.85 -0.35 0.84 0.1608 1.3
## Y 0.86 -0.20 0.78 0.2230 1.1
## Z 0.82  0.56 0.99 0.0063 1.8
## 
##                        PC1  PC2
## SS loadings           2.13 0.48
## Proportion Var        0.71 0.16
## Cumulative Var        0.71 0.87
## Proportion Explained  0.82 0.18
## Cumulative Proportion 0.82 1.00
## 
## Mean item complexity =  1.4
## Test of the hypothesis that 2 components are sufficient.
## 
## The root mean square of the residuals (RMSR) is  0.11 
##  with the empirical chi square  7.66  with prob <  NA 
## 
## Fit based upon off diagonal values = 0.96
summary(fit)
## 
## Factor analysis with Call: principal(r = dat, nfactors = 2, rotate = TRUE)
## 
## Test of the hypothesis that 2 factors are sufficient.
## The degrees of freedom for the model is -2  and the objective function was  0.33 
## The number of observations was  100  with Chi Square =  32  with prob <  NA 
## 
## The root mean square of the residuals (RMSA) is  0.11
biplot(fit)

library(psych)
dat <- d
fit <- principal(dat, nfactors=3, rotate=TRUE)#バリマックス回転
## Specified rotation not found, rotate='none' used
fit
## Principal Components Analysis
## Call: principal(r = dat, nfactors = 3, rotate = TRUE)
## Standardized loadings (pattern matrix) based upon correlation matrix
##    PC1   PC2   PC3 h2       u2 com
## X 0.85 -0.35  0.40  1  0.0e+00 1.8
## Y 0.86 -0.20 -0.47  1 -1.8e-15 1.7
## Z 0.82  0.56  0.08  1 -2.2e-16 1.8
## 
##                        PC1  PC2  PC3
## SS loadings           2.13 0.48 0.39
## Proportion Var        0.71 0.16 0.13
## Cumulative Var        0.71 0.87 1.00
## Proportion Explained  0.71 0.16 0.13
## Cumulative Proportion 0.71 0.87 1.00
## 
## Mean item complexity =  1.8
## Test of the hypothesis that 3 components are sufficient.
## 
## The root mean square of the residuals (RMSR) is  0 
##  with the empirical chi square  0  with prob <  NA 
## 
## Fit based upon off diagonal values = 1
summary(fit)
## 
## Factor analysis with Call: principal(r = dat, nfactors = 3, rotate = TRUE)
## 
## Test of the hypothesis that 3 factors are sufficient.
## The degrees of freedom for the model is -3  and the objective function was  0 
## The number of observations was  100  with Chi Square =  0  with prob <  NA 
## 
## The root mean square of the residuals (RMSA) is  0
biplot(fit)

3 主成分分析の理論.

\(w =\begin{pmatrix} w_{1}\\ \vdots\\ w_{p} \end{pmatrix}~,~~|w|=1\)　、　確率変数 \(X =\begin{pmatrix} X_{1}\\ \vdots\\ X_{p} \end{pmatrix}\)　は平均０に基準化されているとする。

新変量 \(X'w=w_1X_1+\cdots w_pX_p\) とすると.

　　\(E(X'w)=0\)　　\(V(X'w)=E(w'XX'w)=w'E(XX')w=w'\Sigma w\).

wに対して定まる新変量の分散 \(V(X'w)\) は \(w\)の２次形式になっている。
新変量は Xのw方向への正射影である。
新変量の分散が最大となる w を考える.

２次形式の微分公式.
　　　\(\displaystyle{ \frac{\partial}{\partial \mathbf{x}} \mathbf{x}^T\mathbf{A}\mathbf{x} = (\mathbf{A}+\mathbf{A}^T)\mathbf{x} }\)　　　を用いる.

\(|w|=1\) のもとで \(w'\Sigma w\) を最大にするには、

　　　\(L=w'\Sigma w-\lambda(w'w-1)\) を \(w\) で微分する。（ラグランジュの未定乗数法）
二次形式の微分　\(\dfrac{\partial }{\partial w}w'\Sigma w =(\Sigma+\Sigma^T)w= 2 \Sigma w\) を用いると.

　　　\(\dfrac{\partial L}{\partial w} = 2 \Sigma w - 2 \lambda w=0\)　　　\(\therefore \Sigma w =\lambda w\).

よって、

  • \(w\) は \(\Sigma\) の固有ベクトルとして求まる.

  • \(w'\Sigma w=w'\lambda w=\lambda |w|=\lambda\) だから、新変数の分散は固有値 \(\lambda\) に等しい。

固有値と固有ベクトルの求め方については、
当サイト. 2変量の主成分分析 　 固有値と対角化 を参照。

2変量の主成分分析では、

　　\(\Sigma= \begin{pmatrix} s_{11} & s_{12}\\ s_{12}& s_{22} \end{pmatrix}\)　　　\(\begin{pmatrix} s_{11} & s_{12}\\ s_{12}& s_{22} \end{pmatrix} \begin{pmatrix} a_{1}\\ a_{2} \end{pmatrix} = \lambda \begin{pmatrix} a_{1}\\ a_{2} \end{pmatrix}\)　について、

  • 2つの固有値のうち大きいほうを \(\lambda_{(1)}\)、そのときの固有ベクトルを \(\begin{pmatrix} a_{1(1)}\\ a_{2(1)} \end{pmatrix}\)
  • 2つの固有値のうち大きいほうを \(\lambda_{(2)}\)、そのときの固有ベクトルを \(\begin{pmatrix} a_{1(2)}\\ a_{2(2)} \end{pmatrix}\)

とすると.

　　\(\begin{pmatrix} a_{1(1)}&a_{1(2)}\\ a_{2(1)}&a_{2(2)} \end{pmatrix} ^{-1} \begin{pmatrix} s_{11} & s_{12}\\ s_{12}& s_{22} \end{pmatrix} \begin{pmatrix} a_{1(1)}&a_{1(2)}\\ a_{2(1)}&a_{2(2)} \end{pmatrix} = \begin{pmatrix} \lambda_{(1)}&0\\ 0&\lambda_{(2)} \end{pmatrix}\).

直行行列によって対角化できるが、直行行列の性質から.

　　\(\begin{pmatrix} a_{1(1)}&a_{2(1)}\\ a_{1(2)}&a_{2(2)} \end{pmatrix} \begin{pmatrix} s_{11} & s_{12}\\ s_{12}& s_{22} \end{pmatrix} \begin{pmatrix} a_{1(1)}&a_{1(2)}\\ a_{2(1)}&a_{2(2)} \end{pmatrix} = \begin{pmatrix} \lambda_{(1)}&0\\ 0&\lambda_{(2)} \end{pmatrix}\).

対角成分は２次形式の形になっている。