相関がある二変量正規乱数を発生させてみる

　あっ, どーも僕です.

　やたら同級生(女)に相談されるのでチャンスかと勘違いしちゃいました(・∀・)友達一人減ったかしら(・∀・)

相関がある二変量正規乱数

　今日家に帰ってきたら中古で2500円だったのでついポチってしまいまった本が届いていたので, さっそく少し試してみました.
　取り組んだ話題は二変量正規乱数発生です. 本書には線形変換とギブスサンプリングの二通りが掲載されていましたので, 両者をRで組んでみました.

統計的データ解析のための数値計算法入門 (統計ライブラリー)

線形変換

　一変量標準正規分布に従う乱数は発生させられるとして, 互いに独立な確率変数　　を要素に持つ確率ベクトル　　を　
　
となる行列

によって変換します. つまり,

によって目的の乱数が得られます. Rで実装して10万個発生させてみた結果とコードがこちらです. 相関係数は0.5, 0.9です. ちょっと見にくいですが, corが得られた乱数の相関係数です.

# ************************************************
# 相関係数ρをもつ二変量標準正規分布乱数の発生
# ************************************************
# その１
# 標準正規分布乱数を線形変換
r_st_bibnorm <- function (N = 100, rho = .5) 
{
  Z <- matrix (c (rnorm (N), rnorm (N)), nrow = N, ncol = 2)
  A <- matrix (c (1, 0,
                  rho, sqrt (1 - rho ** 2)), nrow = 2, ncol = 2, byrow = TRUE)
  X <- t( tcrossprod (A, Z))
  return (list (values = X, cor = cor (X[ , 1], X[ , 2])))
}

ギブスサンプリング

　ギブスサンプリングによっても乱数が得られます. 定義とかは飛ばしますが, ご存知のように一般の多変量正規分布は次のような形になっています.

　この式のパラメータを

　になるようにしてちょこちょこっろ式変形すると, 片方の確率変数が従う条件付き確率分布は

　になります. やってみればわかりますが, もう片方も同じ形になります. つまり, ギブスサンプリングは次の繰り返しになります.

Rで実装して10万個発生させてみた結果とコードがこちらです. 相関係数は0.5, 0.9です. ちょっと見にくいですが, corが得られた乱数の相関係数です.

# その2
# ギブスサンプリング
gibr_st_bibnorm <- function (N = 120, burn.in = 20, rho = .5)
{
  ini_value <- rnorm (1)
  rtn <- matrix(0, nrow = N, ncol = 2)
  rtn[1, 1] <- ini_value
  SD = sqrt (1 - rho ** 2)
  
  for (i in seq_len(N)[-1]) {
    rtn[i, 2] <- rnorm (1, mean = rho * rtn[i-1, 1], sd = SD)
    rtn[i, 1] <- rnorm (1, mean = rho * rtn[i, 2], sd = SD)
  }
  rtn <- rtn[-(seq_len (burn.in)), ]
  return ( list (values = rtn, cor = cor(rtn[ , 1], rtn[ , 2])))
}

線形変換による一般形

　本書には記載されていませんでしたが, 標準正規分布でなくて一般の正規分布による場合も線形変換で簡単に発生させることができます. 前述の例で既にうまくできているので, 改めて示す必要もないですが, 参考までにRで実装して10万個発生させてみた結果とコードがこちらです. 相関係数は 0.9, -0.9です. ちょっと見にくいですが, corが得られた乱数の相関係数です.
　標準偏差や平均を少しずらしています.

# ************************************************
# その3
# 標準正規分布乱数を線形変換
r_bibnorm <- function (N = 100, rho = .5, 
                       mu1 = 0, mu2 = 0, 
                       sigma1 = 1, sigma2 = 1) 
{
  Z <- matrix (c (rnorm (N), rnorm (N)), nrow = N, ncol = 2)
  A <- matrix (c (sigma1 , 0,
                  rho * sigma2, sqrt (1 - rho ** 2)) * sigma2, 
               nrow = 2, ncol = 2, byrow = TRUE)
  X <- t( tcrossprod (A, Z) + c(mu1, mu2)) 
  return (list (values = X, cor = cor (X[ , 1], X[ , 2]), sd1 = sd(X[,1]), sd2 = sd(X[,2])))
}