多変量正規分布

R 備忘録

論文でm次元ガウス分布を例にして評価実験が行われていたので、それを試してみるために色々と調べていたら、やっぱりありました。自分でゴリゴリ書くのが若干めんどくさかったのでありがたいですね。

多変量解析用パッケージのMASSのmvrnormを使います。

ポイントは引数

  • 第１引数に作りたいデータの個数
  • 第２引数に各次元の平均（ベクトル）
  • 第３引数に共分散行列（行列）

を与えてやります。

第３引数が若干混乱しそうですが、共分散行列は対角成分に分散が、それ以外は例えば1行2列目には1番目と2番目の次元の相関関係が入ります。つまり、各次元が独立している状況を考えたければ単位行列を与えればよいわけです。
第２引数の長さと、共分散行列なので当然正方行列である第３引数の行と列は等しくないといけません。

論文では、N(0,I)とあったので、第２引数は全て０のベクトル、第３引数は第２引数の長さと同じ行と列を持つ単位行列を用意します。

コード

前置きが長くなりましたが

# 多次元ガウス分布に従う乱数を発生させる
library(MASS)

NUM <- 10
DIM <- 5

# 平均0, 分散共分散行列Iとする
mu <- rep(0, DIM)
Sigma <- diag(DIM)

mgauss <- mvrnorm(NUM, mu, Sigma)

これで5次元のガウス分布に従うデータができます。

> mgauss
            [,1]        [,2]       [,3]       [,4]       [,5]
 [1,] -0.5182721  1.32180396  1.1126628  2.4263074  1.0394240
 [2,] -1.8095158 -0.77671453 -1.4647852 -0.3489059 -0.6665914
 [3,]  1.7219589  0.43644475  0.8682914 -0.6138239  1.0527091
 [4,]  0.6804582  2.37529984 -1.2477997 -1.0305826 -0.6582892
 [5,] -2.5596450 -0.03711889 -0.5432647 -0.9534160 -1.2563976
 [6,] -0.2712489 -0.03383668  1.4660665 -0.2652364  1.3676385
 [7,]  0.9598914 -2.38425174  2.2059020 -0.4894067 -2.0569321
 [8,]  0.6628937 -2.23163867  1.7417454  1.0441785  0.7401633
 [9,] -2.0744806  2.60568054  0.6914947  2.1207627 -1.7053113
[10,]  0.7704964 -1.05948353  0.7825592 -0.6031393 -0.6248058

こんな感じの結果になりました。
これが正しいのか調べてみます。

> apply(mgauss, 2, ave)[1,]
[1] -0.24374639  0.02161850  0.56128724  0.12867379 -0.27683925
> apply(mgauss, 2, sd)
[1] 1.462886 1.724264 1.246277 1.267219 1.239451

という感じで、全ての次元がだいたい平均0、標準偏差1くらいになっていますね。

当然変数間の相関関係も気になるので相関係数を求めるcor()を使って調べてみます。

> for (i in 1:(DIM-1)) {
+   for (j in (i+1):DIM) {
+     print(paste(i,j,cor(mgauss[,i], mgauss[,j])))
+   }
+ }
[1] "1 2 -0.297097215983719"
[1] "1 3 0.458390597332640"
[1] "1 4 -0.251747041427307"
[1] "1 5 0.353885520049160"
[1] "2 3 -0.429207042672545"
[1] "2 4 0.293945489581062"
[1] "2 5 0.0165658231607912"
[1] "3 4 0.359103216193921"
[1] "3 5 0.229710900663509"
[1] "4 5 0.178124336924986"

ということで若干相関関係出とるやないかい！

となるのですが、よく考えると今データは10個と統計量を求めるにはちょっと少ないですね。
データ数を1000にしたときの各次元間の相関係数です。

[1] "1 2 0.00150376288361690"
[1] "1 3 -0.0457721539976883"
[1] "1 4 -0.00970859916782877"
[1] "1 5 0.00728152125143677"
[1] "2 3 -0.0115930008821253"
[1] "2 4 0.00541264044206026"
[1] "2 5 0.0231994762339870"
[1] "3 4 -0.00260412182207902"
[1] "3 5 -0.0138434483251586"
[1] "4 5 -0.0500139428051309"

とほぼ0になっており、次元間が無相関であるという感じになってますね。