こんにちはたくまろです。今回は自分の専門分野である、高次元統計解析について紹介したいと思います。データサイエンスやビックデータの解析などで耳にしたことがあるかもしれません。とは言ってもまだまだ新しい分野です。
ここでは、高次元データの主成分分析についてピックアップしたいと思います。特にここでは、ノイズ掃き出し法（Noise-Reduction Methodology;NRM）について紹介します。これは名前の通り、高次元におけるノイズを除去する方法です。また、この記事の内容は僕のバイブルである「高次元の統計学」に基づいています。もっと詳しく知りたい方は、読んでみてください！
 |
高次元の統計学 （統計学One Point 11） [ 青嶋 誠 ] 価格:2,420円 |
はじめに高次元データについて紹介している
【高次元統計解析入門】高次元データのメリット・デメリット【データサイエンス】
を読んでくれると、分かりやすいかもしれません。
主成分分析の問題点について
高次元主成分分析の問題点については下記の記事でまとめています。ご覧になってない方は、先にこちらをよむと何やっているか分かると思います。
ノイズ掃き出し法(NRM)
ノイズ掃き出し（Noise-Reductiob Methodology;NRM）はノイズ項である\( \frac{\sum_{j=i+1}^{d} \lambda_j}{(n-1)\lambda_i}\)を推定して、差し引く方法です。具体的には、次のような式で固有値\(\tilde{\lambda}_j\)を推定します。
[box class=”box5″]
$$\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)$$
[/box]
\(\tilde{\lambda}_j\)は\(\hat{\lambda}_j\)から、ノイズ項を推定したものを差し引いてできた推定量です。ここから固有ベクトルを次のように推定します。
[box class=”box5″]
$$\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)$$
[/box]
ただし、\(\hat{\boldsymbol{h}}_j\)は標本共分散行列\(\boldsymbol{S}_n\)の固有ベクトルで、\(\boldsymbol{u}_j\)は双対標本共分散行列\(\boldsymbol{S}_D\)の固有ベクトルです。
ノイズ掃き出し法(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の方が比の値が次元数が高くなるにつれて１に近づいているのがわかりますね。
コメント