2019-07-16

Rで見る主成分分析と特異値分解の関係

多次元主成分分析を行うにあたって、主成分分析と特異値分解の関係を復習してみた。

まず、データを作成する。

# x：1000個の正規分布N(0, 0.2)に従う乱数
# y : xを３倍し、そこにN(0, 0.3)に従うノイズを加えたもの

# x2, y2: x, yそれぞれについて平均を引いて標準偏差で割って標準化したもの

# x3, y3: x2, y2を大きさ１になるように規格化

# z2: cbind(x3, y3)で作成した1000 x 2の行列

x <- rnorm(1000, 0, 0.2)
y <- 3*x + rnorm(1000, 0, 0.3)
mx <- mean(x)
sdx <- sd(x)
my <- mean(y)
sdy <- sd(y)
x2 <- (x - mx)/sdx
y2 <- (y - my)/sdy

x2 %*% x2
[,1]
[1,] 999

y2 %*% y2
[,1]
[1,] 999

x3 <- x2/sqrt(999)
y3 <- y2/sqrt(999)

z2 <- cbind(x3, y3)

plot(x3, y3)

で確認できるように、第一主成分だけでほぼ説明できるデータである。

z2を特異値分解する。

# w2: z2のSVD
# w2$u 1000 x 2 主成分得点行列

w2 <- svd(z2)

w2の中身を見る。

str(w2)
List of 3
$ d: num [1:2] 1.376 0.327
$ u: num [1:1000, 1:2] 0.00732 0.01653 -0.00181 0.03669 -0.02696 ...
$ v: num [1:2, 1:2] 0.707 0.707 0.707 -0.707

# 主成分と元データの相関係数
cor(w2$u[,1], x3)
[1] 0.972893
cor(w2$u[,1], y3)
[1] 0.972893
cor(w2$u[,2], x3)
[1] 0.2312559
cor(w2$u[,2], y3)
[1] -0.2312559

# w2$uの直交性を利用して、データからw2$uを消して、因子負荷行列を作成
t(w2$u) %*% z2
x3 y3
[1,] 0.9728930 0.9728930
[2,] 0.2312559 -0.2312559

# 特異値の対角行列diad(w2$d)とt(w2$v)から、因子負荷行列を作成
diag(w2$d) %*% t(w2$v)
[,1] [,2]
[1,] 0.9728930 0.9728930
[2,] 0.2312559 -0.2312559

diag(w2$d) %*% t(w2$v) 2 x 2 因子負荷行列
因子負荷は主成分と元データの相関であることが確認できる。

主成分軸の表現
第一主成分のプロット
a <- seq(-0.2,0.2, 0.001)
b1 <- a * w2$v[1,1]
plot(x3, y3, xlim=c(-0.2, 0.2), ylim=c(-0.1, 0.15))
par(new=T)
plot(b1, b2, xlim=c(-0.2, 0.2), ylim=c(-0.1, 0.15), ty='l', col=rgb(1,0,0))b2 <- a * w2$v[2,1]

主成分分析では、w2＄vが相関係数行列の固有値問題を解いた場合の
固有ベクトルとして得られる。
vを主成分得点とした時の負荷行列を計算　
これが元データを回転させた主成分軸上の座標を与える
e <- w2$u %*% diag(w2$d)
max(e[,1])
[1] 0.174004
min(e[,1])
[1] -0.1284722
max(e[,2])
[1] 0.03881079
min(e[,2])
[1] -0.02979512
plot(e[,1],e[,2],xlim=c(-0.13,0.18), ylim=c(-0.13,0.18))

直交性を利用して、t(w2$v)を打ち消して、因子負荷行列を求める
上と同じプロットが得られる。
d <- z2 %*% w2$v
max(d[,1])
[1] 0.174004

max(d[,2])
[1] 0.03881079
min(d[,1])
[1] -0.1284722
min(d[,2])
[1] -0.02979512
plot(d[,1],d[,2],xlim=c(-0.13,0.18), ylim=c(-0.13,0.18))

因子負荷行列が主成分軸上での座標を表している。

2019-07-16

jagsで潜在混合モデル

ベイス統計で実戦モデリング（北大路書房）の第6章　潜在混合モデルの練習で、次のコードを作成してみた。

８個のデータが４個づつ違う正規分布に従うとして、データがどちらの正規分布に従うか、またその時の正規分布の平均を推測する。ただし、標準偏差はどちらも同じであるとする。

モデルは次のように作成した。

model1.txt

model {
for (i in 1: 8) {
z[i] ~ dbern(0.5)
}
lambda ~ dgamma(0.001, 0.001)
psi ~ dbeta(1,1)
m1 <- psi*0.1 + 0.48
phi ~ dbeta(1,1)
m2 <- phi*0.1 + 0.48
for (i in 1:8) {
mx[i] <- equals(z[i],0)*m1 + equals(z[i],1)*m2
k[i] ~ dnorm(mx[i], lambda)
}
}

psiとphiが正規分布の平均値、lambdaが標準偏差

zには１あるいは０が入り、どちらのグループかをこれで区別する。

以下がRのコードである。

rm(list=ls())

setwd("directory指定")

library(R2jags)

k <- c(30,26,11,45,28,10,9,8)

data <- list("k")

x <- c(rep(0,5000),rep(1,5000))

myinits <- list(list(z=sample(x, 8), psi <- 0.5, phi <- 0.5, lambda <- 0.1))

parameters <- c("z", "m1", "m2")

samples <- jags(data, inits=myinits, parameters,model.file="model1.txt", n.chains=1, n.burnin=1,n.iter=10000, n.thin=1, DIC=T)

2019-06-18

rTensorで多次元主成分分析をやってみた

RのパッケージのrTensorを使って多次元主成分分析を行ってみた。

多次元主成分分析のやり方は、

http://www.iri.pref.kumamoto.jp/sgk/2011/cdrom/thesis/04_hall/04_session/444.pdf

からダウンロードできる日本語文献のpdfファイル　”多次元主成分分析による医療データの特徴分析”を参考にした。ただし、記述に間違いがあるようなので、同じ著者らによる英語論文"Application of Multi-Dimensional Principal Component Analysos to Medical Data"

https://pdfs.semanticscholar.org/280e/b22f667a781b5c6b19ddc8e270db111daf6b.pdf

も参考にした。具体的には、日本語論文の(2)式の二階テンソルの平均の計算の説明にある総和をnでわるという処理の意味がわからなかったのだが、英語論文の記述から、日本語文献の記述は間違いで、何について平均をとれば良いのかがわかった（英語論文　式(1))。また、日本語文献では中心化処理しかされていないが、英語論文では標準偏差で割って正規化されていた(英語論文　式(2),(3))。また、日本語文献の式(4), (5)の処理の意味が理解できなかったのだが、英語論文　式(4), (5)から、日本語文献の記述は誤りだとわかった。

HOSVDの処理や意味については

HOSVDとHOOIによる画像の低ランク近似 - Qiita

HOSVDによる画像の低ランク近似ステップ・バイ・ステップ - Qiita

がわかりやすかった。日本語文献には、各モードの主成分の寄与率の計算式は書かれているのだが（式(7)）、特異値をどこから得るのかがわからなかった。HOSVDの処理では、テンソルを行列に変換してSVDを適用しているので、ここから得ているのだろうと推測された。rTensorのhosvdの出力には、特異値の情報は含まれていなさそうなので、寄与率を計算するには、

HOSVDによる画像の低ランク近似ステップ・バイ・ステップ - Qiita

のpythonのコードをもとに自作する必要がありそう。rTensorのコードはGitHubで公開されているようなので、hosvdを書きなして特異値を出力できる関数を作成しても良いかもしれない。

https://github.com/cran/rTensor/blob/master/R/rTensor_Decomp.R