多次元正規分布を作る

# 多次元正規分布乱数を作る関数を作っておく
library(MASS)
library(lqmm)
# 適当に対称行列を作って正定値行列を作る関数も作っておく
# 多次元正規分布は中心座標と分散共分散行列とで決まる
# 分散共分散行列のために正定値行列が必要
my.make.positive.definite <- function(n.v){# n.vは変数の個数
  sigma <- matrix(rnorm(n.v^2),n.v,n.v)
  sigma <- sigma + t(sigma)
  sigma <- make.positive.definite(sigma)
  return(sigma)
}
# n.sampleは発生乱点数、n.vは変数の個数(次元)
# mは中心座標、sigmaはn.v x n.v 正定値行列
my.mvrnorm <- function(n.sample,n.v,m=NULL,sigma=NULL){
  if(is.null(m)){
    m <- rnorm(n.v)
  }
  if(is.null(sigma)){
    sigma <- my.make.positive.definite(n.v)
  }
  return(mvrnorm(n.sample,m,sigma))
}


# サンプル数
n.s <- 1000

# 説明変数
## 説明変数の数
n.v <- 5
## 多次元正規分布とする
M <- rnorm(n.v) # 重心
S <- my.make.positive.definite(n.v) # Sigma行列(分散共分散行列のようなもの)
X <- my.mvrnorm(n.s,n.v,M,S)
plot(as.data.frame(X))

library(rgl)


tmp.X <- X[,1:3]
tmp.X <- rbind(tmp.X,rep(max(tmp.X),3),rep(min(tmp.X),3)) # ３軸スケールをそろえるため

plot3d(tmp.X)

# sigma行列を適当に与えて、多次元正規分布(ただしゆがんでいる)を作った

# 点の数をだんだん増やす
n.ss <- c(10,25,50,75,100,250,500,750,1000,5000,10000)

diff.M <- rep(0,length(n.ss))
for(i in 1:length(n.ss)){
	n.s <- n.ss[i]
	X <- my.mvrnorm(n.s,n.v,M,S)

	# この乱点のラグビーボール状態を調べるのは、PCA
	# 標本重心
	M.hat <- apply(X,2,mean)
	# 真値と比べる
	M.hat
	M

	pca.out <- prcomp(X)
	# 長軸からsdが並んで出る
	# sdを対角成分にした行列を作り

	K <- diag(pca.out$sdev)
	# 回転行列も出ているから、それで回す
	R <- pca.out$rotation

	P <- R %*% K

	Q <- P %*% t(P) # var-covとは各軸の内積のようなもの
	S

	diff.M[i] <- sum((Q-S)^2)

}

plot(log(n.ss),diff.M,type="l")