【高次元統計解析入門】高次元データの主成分分析(PCA)【ノイズ掃き出し法(NRM)】


こんにちはたくまろです。今回は自分の専門分野である、高次元統計解析について紹介したいと思います。データサイエンスビックデータの解析などで耳にしたことがあるかもしれません。とは言ってもまだまだ新しい分野です。

ここでは、高次元データの主成分分析についてピックアップしたいと思います。特にここでは、ノイズ掃き出し法(Noise-Reduction Methodology;NRM)について紹介します。これは名前の通り、高次元におけるノイズを除去する方法です。また、この記事の内容は僕のバイブルである「高次元の統計学」に基づいています。もっと詳しく知りたい方は、読んでみてください!

高次元の統計学 (統計学One Point 11) [ 青嶋 誠 ]

価格:2,420円
(2021/3/16 09:50時点)
感想(0件)

はじめに高次元データについて紹介している

【高次元統計解析入門】高次元データのメリット・デメリット【データサイエンス】

を読んでくれると、分かりやすいかもしれません。

主成分分析の問題点について

高次元主成分分析の問題点については下記の記事でまとめています。ご覧になってない方は、先にこちらをよむと何やっているか分かると思います。

【高次元統計解析入門】高次元データの主成分分析(PCA)【主成分分析の問題点】

ノイズ掃き出し法(NRM)

ノイズ掃き出し(Noise-Reductiob Methodology;NRM)はノイズ項である\( \frac{\sum_{j=i+1}^{d} \lambda_j}{(n-1)\lambda_i}\)を推定して、差し引く方法です。具体的には、次のような式で固有値\(\tilde{\lambda}_j\)を推定します。

$$\tilde{\lambda}_j=\hat{\lambda}_j – \frac{ {\rm tr}(\boldsymbol{S}_D)-\sum_{i=1}^{j} \hat{\lambda}_i}{n-1-j}\ \ (j=1,..,n-2)$$

\(\tilde{\lambda}_j\)は\(\hat{\lambda}_j\)から、ノイズ項を推定したものを差し引いてできた推定量です。ここから固有ベクトルを次のように推定します。

$$\tilde{\boldsymbol{h}}_j = \frac{(\boldsymbol{X}-\bar{\boldsymbol{X}})\hat{\boldsymbol{u}}_j}{\sqrt{(n-1)\tilde{\lambda}_j}}= \sqrt{\frac{\hat{\lambda}_j}{\tilde{\lambda}_j}}\hat{\boldsymbol{h}}_j\ \ (j=1,…,n-1)$$

ただし、\(\hat{\boldsymbol{h}}_j\)は標本共分散行列\(\boldsymbol{S}_n\)の固有ベクトルで、\(\boldsymbol{u}_j\)は双対標本共分散行列\(\boldsymbol{S}_D\)の固有ベクトルです。

IMAGE

イメージはこんな感じです。\(\sqrt{\frac{\hat{\lambda}_i}{\tilde{\lambda}_i}}>1\)なので、\(\tilde{\boldsymbol{h}}_i\)のノルムは\(1\)以上です。もともと大きさ\(1\)だった、\(\hat{\boldsymbol{h}}_j\)()をすこし長くして\(\boldsymbol{h}_j\)(黒)との内積が1になる調整したものが\(\tilde{\boldsymbol{h}}_i\)()です。

ノイズ掃き出し法(NRM)のRコード

RコードによってNRMを表現すると以下のようになります。ただし入力はデータ行列\(X\)です。

NRM <- function(X){
d <- dim(X)[1]
n <- dim(X)[2]
r <- min(n-2, d)
Mean <- apply(X, 1, mean)
X <- sweep(X, 1, Mean, '-')
Sd <- t(X) %*% X / (n - 1)
eig <- eigen(Sd)
dualval <- eig$values[1:r]
dualvec <- eig$vectors
nrmval <- numeric(r)
nrmvec <- matrix(0, d, r)
nrmscore <- matrix(0, n, r)
for (i in 1:r){
nrmval[i] <- dualval[i] - (sum(diag(Sd)) - sum(dualval[1:i])) / (n - i - 1)
nrmvec[, i] <- X %*% dualvec[, i] / sqrt((n - 1) * nrmval[i])
for (j in 1:n){
nrmscore[j, i] <- dualvec[j, i] * sqrt(n * nrmval[i])
}
}
return(list(values=nrmval, vectors=nrmvec, scores=nrmscore))
}

高次元主成分分析の比較

共分散行列の固有値の推定を\(\hat{\lambda}_j\)と\(\tilde{\lambda}_j\)で比較します。\(\hat{\lambda}_j/{\lambda}_j\)と\(\tilde{\lambda}_j/{\lambda}_j\)の比が1になるとき正しく推定出来ていることになります。\(\lambda_1=d^{2/3},\ \lambda_2=d^{1/2}\)として推定します。

datagene_norm <- function(mu,Sigma,N){
d <- dim(Sigma)[1]
Z <- matrix(0,nrow=d,ncol=N)

for (i in 1:N) {
Z[,i] <- rnorm(d,0,1)
}
gamma<- matrix(0,d,d)
diag(gamma) <- eigen(Sigma)$value^(1/2)
H <- eigen(Sigma)$vector

return(H%*%gamma%*%Z+mu)
}

L <- 1000
M <- 7
val11 <- matrix(0,M,L)
val21 <- matrix(0,M,L)
val12 <- matrix(0,M,L)
val22 <- matrix(0,M,L)

for (i in 3:10) {
d <- 2^i
N <- 2*ceiling(d^(1/2))
Sigma <- diag(d)
Sigma[1,1] <- d^(2/3)
Sigma[2,2] <- d^(1/2)
mu <- numeric(d)
X_ <- datagene_norm(mu,Sigma,N*L)
for (l in 1:1000) {
X <- X_[,(1+N*(l-1)):(N*l)]
X <- X-rowMeans(X)
S <- t(X)%*%X/(N-1)
Model1 <- eigen(S)
Model2 <- NRM(X)
val11[i-2,l] <- Model1$values[1]/d^(2/3)
val21[i-2,l] <- Model2$values[1]/d^(2/3)
val12[i-2,l] <- Model1$values[2]/d^(1/2)
val22[i-2,l] <- Model2$values[2]/d^(1/2)
}
}

l1 <- rep(1,7)
cols <- c("blue", "red", "black")
ltys <- c(1, 4, 1)
pchs <- c(1, 1, 1)

data <- data.frame(rowMeans(val11),rowMeans(val21),l1)
plot(0, 0, type = "n",xlim=c(1,7),ylim=c(0.8,2.0), xlab = "dimension" , ylab = "rate of first eigenvalue" , xaxt="n")
for (i in 1:3) {
if(i<3)points(data[, i], col = cols[i])
lines(data[, i], col = cols[i],lty = ltys[i])
}
name <- c(4:11)
axis(1,1:8,labels=name)

data <- data.frame(rowMeans(val12),rowMeans(val22),l1)
plot(0, 0, type = "n",xlim=c(1,7),ylim=c(0.8,2.0), xlab = "dimension" , ylab = "rate of second eigenvalue" , xaxt="n")
for (i in 1:3) {
if(i<3)points(data[, i], col = cols[i])
lines(data[, i], col = cols[i],lty = ltys[i])
}
name <- c(4:11)
axis(1,1:8,labels=name)

\(d=2^s,\ n=2*d^{1/2}\ (s=4,..,10)\) としています。横軸が\(s\)で縦軸が推定値と真の値との比です。青色が従来の推定量で、赤色がNRMによる推定量です。NRMの方が比の値が次元数が高くなるにつれて1に近づいているのがわかりますね。

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です