【Rコード】多変量正規分布のデータ発生【パッケージなし】

こんにちはたくまろです。今回は僕がRでよく使う関数である、多変量正規分布のデータ発生について紹介したいと思います。今回はRコードで多変量正規分布のデータ発生をパッケージなしで書いています。一般に良く知られている関数

mvrnorm(N,mu,Sigma)

でも出来ますが、次元数が高いと計算速度が遅いので、今回は自分でコードを書きました。標準正規分布を発生させてから、データを変換しているためより早く多変量正規分布のデータ発生させています。高次元の正規分布からデータを速く発生させたい方はご参照ください。

多変量正規分布のデータ発生

\( \boldsymbol{x}_1,…,\boldsymbol{x}_N\sim N_d(\boldsymbol{\mu},\boldsymbol{\Sigma}) \)をRコードで発生させると次の通りです。

###input###
# mu     mean vector
# Sigma  covariance matrix
# N      sample size
###output###
# N*d data matrix (N:sample, d:dimension)

datagene_norm <- function(mu,Sigma,N){
d <- dim(Sigma)[1]
Z <- matrix(0,d,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)
}

引数は平均muと共分散行列Sigmaとサンプル数Nです。

多変量正規分布のデータ発生関数を使ってみる

では、\( \boldsymbol{x}_1,…,\boldsymbol{x}_N\sim N_{100}(\boldsymbol{\mu},\boldsymbol{\Sigma}) \)を実際に発生させてみましょう。

d <- 100 #dimension 
N <- 10 #sample size 
mu <- rep(1,d) #mean vector 
Sigma <- matrix(0,d,d) #covariance matrix 
for(i in 1:d){ for(j in 1:d){ Sigma[i,j] <- 0.7^abs(i-j)  } }

X <- datagene_norm(mu,Sigma,N) #Generat d*N data.
dim(X)
[1] 100 10

という風に出力され、正しい大きさになっています。ただし\(\boldsymbol{X}\)は次元×サンプル数のサイズになっています。

コメントを残す

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