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

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

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)