はじめに

　機械学習や統計学で登場する各種の確率分布について、「計算式の導出・計算のスクラッチ実装・計算過程や結果の可視化」などの「数式・プログラム・図」を用いた解説により、様々な角度から理解を目指すシリーズです。

　この記事では、ガウス分布の乱数についてPythonを使って確認します。

【前の内容】

【他の内容】

【今回の内容】

多次元ガウス分布の乱数生成

　多次元ガウス分布(multivariate Gaussian distribution・多変量正規分布・multivariate Normal distribution)に従う乱数(random number)を生成して、グラフやアニメーションで確認します。この記事では、PythonのMatplotlibライブラリを利用して作図します。
　ガウス分布については「多次元ガウス分布の定義式 - からっぽのしょこ」、グラフ作成については「【Python】2次元ガウス分布の作図 - からっぽのしょこ」、Rを利用する場合は「【R】多次元ガウス分布の乱数生成 - からっぽのしょこ」を参照してください。

　利用するライブラリを読み込みます。

# ライブラリを読込
import numpy as np
from scipy.stats import multivariate_normal
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import matplotlib.colors as colors

サンプリング

　まずは、多次元ガウス分布に従う乱数を作成します。
　この例では、2次元のグラフで可視化するため、次元数を $D = 2$ とします。乱数生成の処理自体は次元数に関わらず行えます。

パラメータの設定

　 $D$ 次元ガウス分布のパラメータ $\boldsymbol{\mu}, \boldsymbol{\Sigma}$ を指定します。

# 次元数を設定
D = 2

# 平均パラメータを指定
mu_d = np.array(
    [0.0, 10.0]
)
print(mu_d)

# 分散パラメータを指定
sigma2_dd = np.array(
    [[9.0, 2.5], 
     [2.5, 4.0]]
)
print(sigma2_dd)

[ 0. 10.]
[[9.  2.5]
 [2.5 4. ]]

　平均ベクトル(実数ベクトル) $\boldsymbol{\mu}$ 、分散共分散行列(正定値行列) $\boldsymbol{\Sigma}$ を指定します。

$\displaystyle \boldsymbol{\mu} = \begin{pmatrix} \mu_1 & \mu_2 \end{pmatrix} ,\ \boldsymbol{\Sigma} = \begin{pmatrix} \sigma_1^2 & \sigma_{1,2} \\ \sigma_{2,1} & \sigma_2^2 \end{pmatrix}$

　各軸の平均(実数) $\mu_d$ 、分散(正の値) $\sigma_d^2 \gt 0$ 、2軸の共分散(実数) $\sigma_{i,j}$ に対応します。

　続いて、設定した値に従う乱数を生成して、グラフを作成していきます。

乱数の生成

　サンプルサイズ $N$ を指定して、ガウス分布の乱数 $\mathbf{x}_n$ を生成します。

# サンプルサイズを指定
N = 10000

# ガウス分布の乱数を生成
x_nd = np.random.multivariate_normal(mean=mu_d, cov=sigma2_dd, size=N)
print(x_nd[:5])

[[ 1.68025021  7.14572129]
 [-0.06327498  8.54718483]
 [ 1.77152552 10.48677074]
 [-2.30829669  7.83406996]
 [ 2.84486092  9.93536888]]

　サンプルサイズ(データ数) $N$ を指定します。
　 $N$ 個のサンプル(乱数)の集合 $\mathbf{x} = \{\mathbf{x}_1, \cdots, \mathbf{x}_N\}$ 、 $\mathbf{x}_n = (x_{n,1}, \cdots, x_{n,D})$ を作成します。
　多次元ガウス分布の乱数は、np.random.multivariate_normal() で生成できます。平均ベクトルの引数 mean に $\boldsymbol{\mu}$ 、分散共分散行列の引数 cov に $\boldsymbol{\Sigma}$ 、サンプルサイズの引数 size に $N$ を指定します。
　 $N$ 行 $D$ 列の2次元配列が出力されます。各行が1つのサンプル $\mathbf{x}_n$ に対応します。

　以上で、乱数を作成できました。続いて、乱数を可視化していきます。

乱数の作図

　確率変数 $\mathbf{x}$ の集計範囲を設定します。

# x軸の範囲を設定
u = 5.0
x_0_size = np.abs(x_nd[:, 0]-mu_d[0]).max()
x_1_size = np.abs(x_nd[:, 1]-mu_d[1]).max()
x_0_size = np.ceil(x_0_size /u)*u # u単位で切り上げ
x_1_size = np.ceil(x_1_size /u)*u # u単位で切り上げ
x_0_min  = mu_d[0] - x_0_size
x_0_max  = mu_d[0] + x_0_size
x_1_min  = mu_d[1] - x_1_size
x_1_max  = mu_d[1] + x_1_size
print(x_0_min, x_0_max)
print(x_1_min, x_1_max)

-15.0 15.0
0.0 20.0

　集計や作図に用いるため、階級値の最小値・最大値を作成します。
　この例では、生成したサンプル x_nd の最小値・最大値を使って、範囲を設定しています。

　2次元ガウス分布の乱数のヒートマップを作成します。

# 生成分布のラベルを作成
tmp_mu_str    = ', '.join([f'{mu:.1f}' for mu in mu_d])
tmp_sgm_0_str = ', '.join([f'{sgm:.1f}' for sgm in sigma2_dd[0, :]])
tmp_sgm_1_str = ', '.join([f'{sgm:.1f}' for sgm in sigma2_dd[1, :]])
param_lbl  = f'$N = {N}, '
param_lbl += f'\\mu = ({tmp_mu_str}), '
param_lbl += f'\\Sigma = \\binom{{{tmp_sgm_0_str}}}{{{tmp_sgm_1_str}}}$'

# サンプルの度数を作図
plt.figure(figsize=(9, 6), dpi=100, facecolor='white')
plt.hist2d(
    x=x_nd[:, 0], y=x_nd[:, 1], 
    bins=30
) # 度数
plt.xlabel('$x_1$') # x軸ラベル
plt.ylabel('$x_2$') # y軸ラベル
plt.title(param_lbl, loc='left') # タイトル
plt.suptitle('multivariate Gaussian distribution', fontsize=20) # 図タイトル
plt.colorbar(label='frequency') # カラーバー
plt.grid() # グリッド線
plt.xlim(xmin=x_0_min, xmax=x_0_max) # x軸の表示範囲
plt.ylim(ymin=x_1_min, ymax=x_1_max) # y軸の表示範囲
plt.show()

　未集計データのヒートマップは、plt.hist2d() で描画できます。データの引数 x, y にサンプルデータ、バーの数の引数 bins に(自動で設定される集計範囲における)バーの数を指定します。

　上の例では、階級数(バーの数)を指定して、階級幅や階級値は自動で設定されます。
　次の例では、階級幅や階級値を指定して、任意のバーを作図します。

　階級数から階級幅を作成します。

# 階級数を指定
class_0_num = 31
class_1_num = 21

# 境界数を設定
bin_0_num = class_0_num + 1
bin_1_num = class_1_num + 1

# 階級幅を計算
bin_0_size = (x_0_max - x_0_min) / (class_0_num-1)
bin_1_size = (x_1_max - x_1_min) / (class_1_num-1)
print(bin_0_size)
print(bin_1_size)

1.0
1.0

　階級数 bin_*_num を用いて、階級幅 bin_*_size を計算します。
　階級値の最小値を $x_{\min}$ 、最大値を $x_{\max}$ として、階級幅(1区間のサイズ) $w$ は、階級値の範囲 $x_{\max} - x_{\min}$ を階級数(区間の数) $m$ から1つ減らした数で割った値 $w = \frac{x_{\max} - x_{\min}}{m-1}$ により求まります。

　または、階級幅から階級数を作成します。

# 階級幅を指定
bin_0_size = 1.0
bin_1_size = 1.0

# 階級数を計算
bin_0_num = ((x_0_max - x_0_min) // bin_0_size).astype(dtype='int') + 1
bin_1_num = ((x_1_max - x_1_min) // bin_1_size).astype(dtype='int') + 1
print(bin_0_num)
print(bin_1_num)

31
21

　階級幅 bin_*_size を用いて、階級数 bin_*_num を計算します。
　階級数(区間の数) $m$ は、階級値の範囲 $x_{\max} - x_{\min}$ を階級幅(1区間のサイズ) $w$ で割った整数に1つ増やした数 $m = \lfloor \frac{x_{\max} - x_{\min}}{w} \rfloor + 1$ により求まります。
　 $\frac{x_{\max} - x_{\min}}{w}$ が割り切れない(整数にならない)場合は、 $w$ を再設定する必要があります。
　作図の都合上、階級数(バーの数)を整数型にしておきます。

　階級値と境界値を作成します。

# 境界値の範囲を設定
bin_0_min = x_0_min - 0.5*bin_0_size
bin_0_max = x_0_max + 0.5*bin_0_size
bin_1_min = x_1_min - 0.5*bin_1_size
bin_1_max = x_1_max + 0.5*bin_1_size
print(bin_0_min, bin_0_max)
print(bin_1_min, bin_1_max)

# 階級値を作成
center_0_vec = np.arange(start=x_0_min, stop=x_0_max+bin_0_size, step=bin_0_size)
center_1_vec = np.arange(start=x_1_min, stop=x_1_max+bin_1_size, step=bin_1_size)
print(center_0_vec[:5])
print(center_1_vec[:5])
print(center_0_vec.shape)
print(center_1_vec.shape)

# 境界値を作成
bin_0_vec = np.arange(start=bin_0_min, stop=bin_0_max+0.5*bin_0_size, step=bin_0_size) # (x+0.5*wはstop引数が未満のための対策用)
bin_1_vec = np.arange(start=bin_1_min, stop=bin_1_max+0.5*bin_1_size, step=bin_1_size) # (x+0.5*wはstop引数が未満のための対策用)
#bin_0_vec = np.hstack([center_0_vec, center_0_vec[-1]+bin_0_size]) - 0.5*bin_0_size
#bin_1_vec = np.hstack([center_1_vec, center_1_vec[-1]+bin_1_size]) - 0.5*bin_1_size
print(bin_0_vec[:5])
print(bin_1_vec[:5])
print(bin_0_vec.shape)
print(bin_1_vec.shape)

-15.5 15.5
-0.5 20.5
[-15. -14. -13. -12. -11.]
[0. 1. 2. 3. 4.]
(31,)
(21,)
[-15.5 -14.5 -13.5 -12.5 -11.5]
[-0.5  0.5  1.5  2.5  3.5]
(32,)
(22,)

　階級値の最小値 $x_{\min}$ から最大値 $x_{\max}$ までを階級幅 $w$ の間隔で分割して、各区間の階級値(境界値の中央値)を作成します。
　境界値(区間全体)の最小値 $x_{\min} - \frac{w}{2}$ から最大値 $x_{\max} + \frac{w}{2}$ までを階級幅 $w$ の間隔で分割して、各区間の境界値(階級値の中央値)を作成します。または、各境界値から階級幅の半分 $\frac{w}{2}$ をスライドした値を作成します。

　階級を指定して、ヒートマップを作成します。

# サンプルの度数を作図
plt.figure(figsize=(9, 6), dpi=100, facecolor='white')
plt.hist2d(
    x=x_nd[:, 0], y=x_nd[:, 1], 
    bins=(class_0_num, class_1_num), 
    range=[(bin_0_min, bin_0_max), (bin_1_min, bin_1_max)]
) # 度数
plt.xlabel('$x_1$') # x軸ラベル
plt.ylabel('$x_2$') # y軸ラベル
plt.title(param_lbl, loc='left') # タイトル
plt.suptitle('multivariate Gaussian distribution', fontsize=20) # 図タイトル
plt.colorbar(label='frequency') # カラーバー
plt.grid() # グリッド線
plt.xlim(xmin=x_0_min, xmax=x_0_max) # x軸の表示範囲
plt.ylim(ymin=x_1_min, ymax=x_1_max) # y軸の表示範囲
plt.show()

　バーの数の引数 bins に階級数、バーの範囲の引数 range に境界値(区間全体)の最小値・最大値を指定します。それぞれ軸ごとに値を設定できます。バーの幅は引数の設定により決まります。

　密度のグラフを作成します。

# サンプルの度数を作図
plt.figure(figsize=(9, 6), dpi=100, facecolor='white')
plt.hist2d(
    x=x_nd[:, 0], y=x_nd[:, 1], 
    bins=(class_0_num, class_1_num), 
    range=[(bin_0_min, bin_0_max), (bin_1_min, bin_1_max)], 
    density=True
) # 密度
plt.xlabel('$x_1$') # x軸ラベル
plt.ylabel('$x_2$') # y軸ラベル
plt.title(param_lbl, loc='left') # タイトル
plt.suptitle('multivariate Gaussian distribution', fontsize=20) # 図タイトル
plt.colorbar(label='density') # カラーバー
plt.grid() # グリッド線
plt.xlim(xmin=x_0_min, xmax=x_0_max) # x軸の表示範囲
plt.ylim(ymin=x_1_min, ymax=x_1_max) # y軸の表示範囲
plt.show()

　density=True を指定すると、密度に変換したヒートマップを描画します。

　ここまでで、作図関数内の集計処理によって、グラフを作成できました。続いては、データフレーム上での集計処理によって、集計結果と共に同様のグラフを作成します。

乱数の集計

　サンプル集合 $\mathbf{x}$ を集計します。

# 度数を集計
smp_freq_grid, bin_0_vec, bin_1_vec = np.histogram2d(
    x=x_nd[:, 0], y=x_nd[:, 1], 
    bins=(class_0_num, class_1_num), 
    range=[(bin_0_min, bin_0_max), (bin_1_min, bin_1_max)]
)
print(bin_0_vec[:5])
print(bin_1_vec[:5])
print(smp_freq_grid[10:15, 10:15])

# 密度を計算
smp_dens_grid = smp_freq_grid / (bin_0_size * bin_1_size * N)
print(smp_dens_grid[10:15, 10:15])

# 階級値を作成
center_0_vec = bin_0_vec[:-1] + 0.5*bin_0_size
center_1_vec = bin_1_vec[:-1] + 0.5*bin_1_size
print(center_0_vec[:5])
print(center_1_vec[:5])
print(len(center_0_vec))
print(len(center_1_vec))

[-15.5 -14.5 -13.5 -12.5 -11.5]
[-0.5  0.5  1.5  2.5  3.5]
[[ 51.  41.  15.   4.   1.]
 [ 96.  62.  34.  14.   2.]
 [146. 108.  61.  19.   5.]
 [225. 167.  99.  28.  10.]
 [278. 202. 130.  51.  12.]]
[[0.0051 0.0041 0.0015 0.0004 0.0001]
 [0.0096 0.0062 0.0034 0.0014 0.0002]
 [0.0146 0.0108 0.0061 0.0019 0.0005]
 [0.0225 0.0167 0.0099 0.0028 0.001 ]
 [0.0278 0.0202 0.013  0.0051 0.0012]]
[-15. -14. -13. -12. -11.]
[0. 1. 2. 3. 4.]
31
21

（ここの説明が抜けていたので後で追記します）

　または、density=True を指定すると、密度を出力します。

# 密度を計算
smp_dens_grid, _, _ = np.histogram2d(
    x=x_nd[:, 0], y=x_nd[:, 1], 
    bins=(class_0_num, class_1_num), 
    range=[(bin_0_min, bin_0_max), (bin_1_min, bin_1_max)], 
    density=True
)
print(bin_0_vec[:5])
print(bin_1_vec[:5])
print(smp_dens_grid[10:15, 10:15])

[-15.5 -14.5 -13.5 -12.5 -11.5]
[-0.5  0.5  1.5  2.5  3.5]
[[0.0051 0.0041 0.0015 0.0004 0.0001]
 [0.0096 0.0062 0.0034 0.0014 0.0002]
 [0.0146 0.0108 0.0061 0.0019 0.0005]
 [0.0225 0.0167 0.0099 0.0028 0.001 ]
 [0.0278 0.0202 0.013  0.0051 0.0012]]

　ガウス分布の乱数のヒートマップを作成します。

# サンプルの度数を作図
plt.figure(figsize=(9, 6), dpi=100, facecolor='white')
plt.pcolormesh(
    center_0_vec, center_1_vec, smp_dens_grid.T, 
    shading='nearest'
) # 度数
plt.xlabel('$x_1$') # x軸ラベル
plt.ylabel('$x_2$') # y軸ラベル
plt.title(param_lbl, loc='left') # タイトル
plt.suptitle('multivariate Gaussian distribution', fontsize=20) # 図タイトル
plt.colorbar(label='frequency') # 密度
plt.grid() # グリッド線
plt.xlim(xmin=x_0_min, xmax=x_0_max) # x軸の表示範囲
plt.ylim(ymin=x_1_min, ymax=x_1_max) # y軸の表示範囲
plt.show()

　集計済みデータのヒートマップは、plt.pcolormesh() で描画できます。第1・第2引数に階級値、第3引数に度数、セルの色付け方の引数 shading に 'nearest' を指定します。

　ガウス分布に従う確率密度を計算します。装飾用の処理です。

# x軸の値を作成
x_0_vec = np.linspace(start=x_0_min, stop=x_0_max, num=251)
x_1_vec = np.linspace(start=x_1_min, stop=x_1_max, num=251)

# 格子点を作成
x_0_grid, x_1_grid = np.meshgrid(x_0_vec, x_1_vec)

# 格子点の座標を整形
x_arr  = np.stack([x_0_grid.flatten(), x_1_grid.flatten()], axis=1)
x_dims = x_0_grid.shape

# ガウス分布の確率密度を計算
gen_dens_grid = multivariate_normal.pdf(
    x=x_arr, mean=mu_d, cov=sigma2_dd
).reshape(x_dims)

　詳しくは「分布の作図」を参照してください。

　密度のグラフを作成します。

# 密度軸の範囲を設定
u = 0.005
dens_min = 0.0
dens_max = max(smp_dens_grid.max(), gen_dens_grid.max())
dens_max = np.ceil(dens_max /u)*u # u単位で切り上げ
print(dens_max)

# サンプルの度数を作図
plt.figure(figsize=(11, 6), dpi=100, facecolor='white')
pcm = plt.pcolormesh(
    center_0_vec, center_1_vec, smp_dens_grid.T, 
    shading='nearest', 
    alpha=0.5, vmin=dens_min, vmax=dens_max, 
) # 密度
cont = plt.contour(
    x_0_grid, x_1_grid, gen_dens_grid, 
    vmin=dens_min, vmax=dens_max, 
    linewidths=1.0, linestyles='--'
) # 確率密度
plt.xlabel('$x_1$') # x軸ラベル
plt.ylabel('$x_2$') # y軸ラベル
plt.title(param_lbl, loc='left') # タイトル
plt.suptitle('multivariate Gaussian distribution', fontsize=20) # 図タイトル
plt.colorbar(mappable=cont, label='density (generator)') # 確率密度
plt.colorbar(mappable=pcm, label='density (random number)') # 密度
plt.grid() # グリッド線
plt.xlim(xmin=x_0_min, xmax=x_0_max) # x軸の表示範囲
plt.ylim(ymin=x_1_min, ymax=x_1_max) # y軸の表示範囲
plt.show()

0.035

　第3引数に密度(ヒストグラム全体に対するバーの面積の割合)を指定します。

　ガウス分布の乱数の3Dヒストグラムを作成します。

# 配色の設定
cmap = cm.viridis # カラーマップ
norm = colors.Normalize(vmin=dens_min, vmax=dens_max) # 正規化

# 格子点の座標を整形
bin_0_grid, bin_1_grid       = np.meshgrid(bin_0_vec[:-1], bin_1_vec[:-1]) # 境界値の左下端
center_0_grid, center_1_grid = np.meshgrid(center_0_vec, center_1_vec)     # 階級値

# サンプルの密度を作図
fig, ax = plt.subplots(
    figsize=(9, 6), dpi=100, facecolor='white', 
    constrained_layout=True, subplot_kw={'projection': '3d'}
)
fig.suptitle('multivariate Gaussian distribution', fontsize=20)
ax.contour(
    X=x_0_grid, Y=x_1_grid, Z=gen_dens_grid, offset=0.0, 
    vmin=dens_min, vmax=dens_max, 
    linewidths=1.0, linestyles=':'
) # 確率密度
ax.plot_wireframe(
    X=x_0_grid, Y=x_1_grid, Z=gen_dens_grid, 
    color='green', linewidth=0.5, linestyle='--', 
    label='generator'
) # 確率密度
ax.bar3d(
    x=bin_0_grid.flatten(), y=bin_1_grid.flatten(), z=0.0, 
    dx=bin_0_size, dy=bin_1_size, dz=smp_dens_grid.T.flatten(), 
    color=[cmap(norm(z)) for z in smp_dens_grid.T.flatten()], alpha=0.3, 
    label='random number'
) # 密度
ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$')
ax.set_zlabel('density')
ax.set_title(param_lbl, loc='left')
ax.legend(prop={'size': 8})
ax.set_zlim(zmin=dens_min, zmax=dens_max)
#ax.view_init(elev=10, azim=-60) # 表示角度
plt.show()