任意の母相関を持つ多変量の生成

例えば、母相関係数が0.9である変数xと変数yを生成したいとする。

google先生で調べると次のようなサイトが次のような資料を見つける。

Rで架空データの発生

なので、言われた通りにやってみる。

rho <- 0.9 
n=1000
z <- rnorm(n,0,1)
e1 <- rnorm(n,0,1)
e2 <- rnorm(n,0,1)

x <- sqrt(rho)*z+sqrt(1-rho)*e1
y <- sqrt(rho)*z+sqrt(1-rho)*e2

で、相関を確認するとこんな感じ。

> cor(x,y)
[1] 0.9003037

散布図を書いて見るとこんな感じ。

しかし、世の中には便利なものがあるもので、これを一発で出来るパッケージがあることを最近知った（正確には教えてもらった）。使うのはMASSというパッケージに入っているmvrnormという関数。最初から入っているらしい。

このサイトに紹介がある。
Colorless Green Ideas:R言語で二変量正規分布を扱う

library(MASS)

平均のベクトルと分散共分散行列を指定すれば、あとは勝手に作ってくれる。

mu <- c(0,0)
sigma <- matrix(c(1,0.9,0.9,1), nrow = 2, ncol = 2)

dat <- mvrnorm(1000, mu = mu, Sigma = sigma)

散布図書いてみると、良さそうな感じですね。

相関係数はこんな感じ。

> cor(dat)
          [,1]      [,2]
[1,] 1.0000000 0.9049674
[2,] 0.9049674 1.0000000

この関数を使えば、2変量に限らずそれぞれの変数間が任意の母相関を持つn変量を簡単に作ることができる。

mu2 <- c(0,0,0)
sigma2 <- matrix(c(1,0.9,0.7,0.9,1,0.6,0.6,0.9,1), nrow = 3, ncol = 3)

dat2 <- mvrnorm(1000, mu = mu2, Sigma = sigma2)

相関行列を見てみましょう。

> cor(dat2)
          [,1]      [,2]      [,3]
[1,] 1.0000000 0.9026331 0.6867450
[2,] 0.9026331 1.0000000 0.5882691
[3,] 0.6867450 0.5882691 1.0000000

大変便利な世の中ですね。