Pythonで多変量正規分布に従う乱数を生成する

ベクトル自己回帰のダミーデータを生成するために、多変量正規分布に従う乱数が必要なので、
Pythonで生成する方法を紹介します。

numpyとscipyにそれぞれ用意されています。
同じ名前の関数だったので、どちらかの実装をもう一方がラップしているのかと思っていたのですが、
引数の微妙な違いなどあり、どうやら個別に実装されているようです。

ドキュメントはそれぞれ次のページにあります。

numpy
numpy.random.multivariate_normal
(この記事の numpy は version 1.16を使っています。 numpy 1.17.0 のリリースノートを見ると、random moduleに変更が加えられており、どうやらこの関数にも影響が出てるようなのでご注意ください。)

scipy
scipy.stats.multivariate_normal

さて、実際に期待値と分散共分散行列を指定してそれぞれ乱数を生成してみましょう。


import numpy as np
from scipy.stats import multivariate_normal

# 期待値と分散共分散行列の準備
mean = np.array([3, 5])
cov = np.array([[4, -1.2], [-1.2, 1]])

# numpy を用いた生成
data_1 = np.random.multivariate_normal(mean, cov, size=200)

# データ型の確認
print(data_1.shape)
# (200, 2)

# 期待値の確認
print(np.mean(data_1, axis=0))
# [3.00496708 4.94669956]

# 分散共分散の確認
print(np.cov(data_1, rowvar=False))
"""
[[ 3.86542859 -1.31389501]
 [-1.31389501  0.93002097]]
"""

# scipyで生成する方法
data_2 = multivariate_normal(mean, cov).rvs(size=200)

# データ型の確認
print(data_2.shape)
# (200, 2)

# 期待値の確認
print(np.mean(data_2, axis=0))
# [2.81459692 5.10444347]

# 分散共分散の確認
print(np.cov(data_2, rowvar=False))
"""
[[ 4.46151626 -1.28084696]
 [-1.28084696  1.06831954]]
"""

それぞれきちんと生成できたようです。

分散共分散行列の正定値性のバリデーションなど細かなオプションを持っていますが、
あまり使う機会はなさそうです。
(きちんと行う場合も、事前に固有値を求めて確認しておけば大丈夫だと思います。)

コメントを残す

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