【高次元統計解析入門】高次元データの主成分分析【クロスデータ行列法（CDM）】

SHARE

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

ここでは、高次元データの主成分分析についてピックアップしたいと思います。特にここでは、クロスデータ行列法（Cross-Data-Matrix Methodology;CDM）について紹介します。また、この記事の内容は僕のバイブルである「高次元の統計学」に基づいています。もっと詳しく知りたい方は、読んでみてください！

高次元の統計学 （統計学One Point　11） [ 青嶋 誠 ]

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

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

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

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

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

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


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

クロスデータ行列法(CDM)

クロスデータ行列法（Cross-Data-Matrix Methodology;CDM）はそもそも、ノイズ項が発生しないような統計量を使った主成分分析の方法です。具体的には次のように統計量を構築します。

  1. \(\boldsymbol{X}=(\boldsymbol{x}_1,…,\boldsymbol{x}_n)\)を2分割し, \(\boldsymbol{X}_1=(\boldsymbol{x}_{11},…,\boldsymbol{x}_{1n_1})\), \(\boldsymbol{X}_2=(\boldsymbol{x}_{21},…,\boldsymbol{x}_{1n_2})\)とおく. ただし, \(n1\)は\(n/2\)以上の最小の整数で\(n_2=n-n_1\)である.
  2. $$\boldsymbol{S}_{D(1)}=\frac{(\boldsymbol{X}_1-\bar{\boldsymbol{X}}_1)^T(\boldsymbol{X}_2-\bar{\boldsymbol{X}}_2)}{\sqrt{(n_1-1)(n_2-1)}},\ \ \boldsymbol{S}_{D(2)}=\boldsymbol{S}_{D(1)}^T$$
    を計算する.
  3. 対称行列\(\boldsymbol{S}_{D(1)}\boldsymbol{S}_{D(2)}\)の固有値\(\acute{\lambda}_1^2,…,\acute{\lambda}_{n_2-1}^2(\leq 0)\)と対応する固有ベクトルを\(\acute{\boldsymbol{u}}_{1s}\)を計算する. 同様に\(\boldsymbol{S}_{D(2)}\boldsymbol{S}_{D(1)}\)の\(\acute{\lambda}_s^2\)に対応する固有ベクトル\(\acute{\boldsymbol{u}}_{2s}\)も計算する.
  4. さらに符号を揃えるために \(\acute{\boldsymbol{u}}_{2s}=sign(\acute{\boldsymbol{u}}_{2s}^T\boldsymbol{S}_{D(2)}\acute{\boldsymbol{u}}_{1s})\acute{\boldsymbol{u}}_{2s}\)として新たに計算する．
  5. $$\acute{\boldsymbol{h}}_{1i} = \frac{(\boldsymbol{X}_1-\bar{\boldsymbol{X}}_1)\acute{\boldsymbol{u}}_{1i}}{\sqrt{(n_1-1)\acute{\lambda}_i}},\ \ \acute{\boldsymbol{h}}_{2i} = \frac{(\boldsymbol{X}_2-\bar{\boldsymbol{X}}_2)\acute{\boldsymbol{u}}_{2i}}{\sqrt{(n_2-1)\acute{\lambda}_i}}\ \ (i=1,…,n_2-1)$$
    として固有ベクトル\(\acute{\boldsymbol{h}}_{i}\)を次で定義する．
    $$\acute{\boldsymbol{h}}_{i}=(\acute{\boldsymbol{h}}_{1i}+\acute{\boldsymbol{h}}_{2i})/||\acute{\boldsymbol{h}}_{1i}+\acute{\boldsymbol{h}}_{2i}||$$

    サンプルに関して2分割することによって、掛け合わせたときにバイアス項を小さくすることが出来ます。

    さらに、クロスデータ行列法を使う場合、分布に正規性を仮定する必要がありません。正規性というのは、正規分布を緩めた条件です。言い換えれば、t分布のような正規性がない分布や、歪んだ分布にも対応させることが出来ます。

    クロスデータ行列（CDM）のRコード

    RコードによってNRMを表現すると以下のようになります。ただし入力はデータ行列\(X\)です。また、データの2分割はさまざまなやり方があるのでrandom=Falseの場合は単純に2分割しています(デフォルト)。random=Trueの場合はランダムに2分割します。

    CDM <- function(X, random='False'){
d <- dim(X)[1]
n <- dim(X)[2]
n1 <- as.integer(ceiling(n/2))
n2 <- n - n1
n12 <- list(n1, n2)
r <- min(n2-1, d)

if (random=='False'){
index <- c(1:n)
Xcdm <- list(X[, 1:n1], X[, (n1+1):n])
} else if (random=='True'){
pi <- matrix(1/n, 1, n)
index <- sample(c(1:n), n, replace=FALSE, pi)
Xcdm <- list(X[, index[1:n1]], X[, index[(n1+1):n]])
}
Mean <- list(apply(Xcdm[[1]], 1, mean), apply(Xcdm[[2]], 1, mean))
for (i in 1:2){
Xcdm[[i]] <- sweep(Xcdm[[i]], 1, Mean[[i]], '-')
}

Sd <- t(Xcdm[[1]]) %*% Xcdm[[2]] / sqrt((n1 - 1) * (n2 - 1))
eig <- svd(Sd)
cdmval <- eig$d[1:r]
crossvec <- list(eig$u, eig$v)

cdmvec <- matrix(0, d, r)
cdmscore <- matrix(0, n, r)

for (i in 1:r){
crossvec[[2]][, i] <- sign(as.numeric(crossvec[[1]][, i] %*% t(Xcdm[[1]]) %*% Xcdm[[2]] %*% crossvec[[2]][, i])) * crossvec[[2]][, i]

#eigenvalue and eigenvector
h <- list(0, 0)
for (j in 1:2){
h[[j]] <- Xcdm[[j]] %*% crossvec[[j]][, i] / sqrt(cdmval[i] * (n12[[j]] - 1))
}

cdmvec[, i] <- as.vector(h[[1]] + h[[2]]) / 2
cdmvec[, i] <- cdmvec[, i] / norm(cdmvec[, i])


#score
r_score <- list(numeric(n1), numeric(n2))
for (j in 1:2){
for (k in 1:n12[[j]]){
r_score[[j]][k] <- crossvec[[j]][k, i] * sqrt(n12[[j]] * cdmval[i])
}
}
for (k in 1:n){
cdmscore[index[k], i] <- c(r_score[[1]], r_score[[2]])[k]
}
}

return(list(values=cdmval, vectors=cdmvec, scores=cdmscore))
}

    高次元主成分分析の性能比較

    t分布を用いて、従来の固有値推定\(\hat{\lambda}_j\)、NRMの固有値推定\(\tilde{\lambda}_j\)、CDMの固有値推定\(\acute{\lambda}_j\)を比較しました。NRM（ノイズ掃き出し法）についてはこちらで詳しく説明しています。

    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))
}

datagene_t <- function(mu,Sigma,df,N){
d <- dim(Sigma)[1]
Sigma <- ((df-2)/df)*Sigma #correct covariance matrix
X <- datagene_norm(rep(0,d),Sigma,N) #generate N samples by N(0,Sigma)
Z <- datagene_norm(rep(0,df),diag(df),N) #generate N sample by N(0,I)
xi <- apply(Z^2,2,sum)/df #generate xi such that conformed Xi square distribution by Z
Y <- X %*% diag(sqrt(1/xi)) #make t by X and xi
t <- Y + mu #move in mu's direction
return(t)
}

L <- 1000
M <- 7

val11 <- matrix(0,M,L)
val21 <- matrix(0,M,L)
val31 <- matrix(0,M,L)

val12 <- matrix(0,M,L)
val22 <- matrix(0,M,L)
val32 <- matrix(0,M,L)


for (i in 4:10) {
d <- 2^i
N <- ceiling(d^(2/3))
Sigma <- diag(d)
Sigma[1,1] <- d^(2/3)
Sigma[2,2] <- d^(1/2)
mu <- numeric(d)
X_ <- datagene_t(mu,Sigma,df=4,N*L)

for (l in 1:L) {
X <- X_[,(1+N*(l-1)):(N*l)]
X <- X-rowMeans(X)
S <- t(X)%*%X/(N-1)
Model1 <- eigen(S)
Model2 <- CDM(X)
Model3 <- NRM(X)
val11[i-3,l] <- Model1$values[1]/d^(2/3)
val21[i-3,l] <- Model2$values[1]/d^(2/3)
val31[i-3,l] <- Model3$values[1]/d^(2/3)

val12[i-3,l] <- Model1$values[2]/d^(1/2)
val22[i-3,l] <- Model2$values[2]/d^(1/2)
val32[i-3,l] <- Model3$values[2]/d^(1/2)
}
}

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

data <- data.frame(rowMeans(val11),rowMeans(val21),l1,rowMeans(val31))
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:4) {
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,rowMeans(val32))
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:4) {
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)

    青：\(\hat{\lambda}_2/{\lambda}_2\)

    紫：\(\tilde{\lambda}_2/{\lambda}_2\)

    赤：\(\acute{\lambda}_2/{\lambda}_2\)

    正規性が仮定できないような状況ではクロスデータ行列法（CDM）が高次元の主成分分析として有効であることが分かりますね。

    コメントを残す

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

    前の記事

    【高次元統計解析入門】高次元データの主成分分析（PCA）【ノイ…

    次の記事

    【簡単】wordpressでTeXの数式を太文字にする方法