読者です 読者をやめる 読者になる 読者になる

多変量正規分布

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

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

ポイントは引数

  • 第1引数に作りたいデータの個数
  • 第2引数に各次元の平均(ベクトル)
  • 第3引数に共分散行列(行列)

を与えてやります。

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

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

コード

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

# 多次元ガウス分布に従う乱数を発生させる
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になっており、次元間が無相関であるという感じになってますね。