Numpyで多変量正規分布を算出＆プロット

numpy では、N次元の正規分布を出力するnumpy.random.normal() がありますが、こちらでは共分散を設定できません。

共分散を指定する場合、N次元の多変量正規分布を 出力できる、numpy.random.multivariate_normal () を使う必要があります。

numpy.random.multivariate_normal(mean, cov, size)

ここで、キーワード引数sizeに出力する配列の shapeを入力します。

サンプルコード

import numpy as np
import matplotlib.pyplot as plt


n_sample = 1000
n_dim = 4

mean = np.random.normal(0.0, 1, n_dim)
cov = np.eye(n_dim)
for k in range(0, n_dim):
    cov[k,k] = np.abs(np.random.normal(0.0, 1))
for n in range(0, n_dim):
    for m in range(0, n):
        r = np.random.normal(0.0, 1)
        cov[m,n] = r
        cov[n,m] = r
x = np.random.multivariate_normal(mean, cov, n_sample).T
print(x.shape)

count = 1
for m in range(0, n_dim):
    for n in range(0, n_dim):
        if m == n:
            plt.subplot(n_dim, n_dim, count)
            hist, _ = np.histogram(x[n], bins=30, range=(-5., 5.))
            plt.bar(np.arange(0, len(hist)), hist, width=1.0, color="m", alpha=0.5)
            plt.xlabel(f"x[{n}]")
        elif m > n:
            plt.subplot(n_dim, n_dim, count)
            plt.scatter(x[m], x[n], s=1, c="c", alpha=0.25)
            plt.xlabel(f"x[{n}]")
            plt.ylabel(f"x[{m}]")
        count += 1

実行結果