はじめに

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

　この記事では、ガウス分布のグラフ作成についてPythonを使って確認します。

【前の内容】

【他の内容】

【今回の内容】

多次元ガウス分布の作図

　多次元ガウス分布(multivariate Gaussian distribution・多変量正規分布・multivariate Normal distribution)のグラフを作成します。この記事では、PythonのMatplotlibライブラリを利用して作図します。
　ガウス分布については「分布の定義式」、確率密度計算については「分布の計算」、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}$

　確率変数 $x_1, x_2$ の作図範囲を設定します。

# x軸の範囲を設定
#x_0_min = -10.0
#x_0_max = 10.0
#x_1_min = 0.0
#x_1_max = 20.0
k = 3.0
u = 5.0
x_0_size  = np.sqrt(sigma2_dd[0, 0]) # 基準値を指定
x_1_size  = np.sqrt(sigma2_dd[1, 1]) # 基準値を指定
x_0_size *= k # 定数倍
x_1_size *= k # 定数倍
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)

# 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)
print(x_0_vec[:5])
print(x_1_vec[:5])

-10.0 10.0
0.0 20.0
[-10.    -9.92  -9.84  -9.76  -9.68]
[0.   0.08 0.16 0.24 0.32]

　範囲を指定して、確率変数がとり得る値(実数) $x_1, x_2$ の数値ベクトルを作成します。
　この例では、指定したパラメータ mu_d, sigma2_dd を使って、平均を中心に標準偏差の倍数 $\mu_d \pm k \sigma_d$ の範囲を設定しています。

　格子点 $(x_1, x_2)$ を作成します。

# 格子点を作成
x_0_grid, x_1_grid = np.meshgrid(x_0_vec, x_1_vec)
print(x_0_grid[:5, :5])
print(x_1_grid[:5, :5])

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

[[-10.    -9.92  -9.84  -9.76  -9.68]
 [-10.    -9.92  -9.84  -9.76  -9.68]
 [-10.    -9.92  -9.84  -9.76  -9.68]
 [-10.    -9.92  -9.84  -9.76  -9.68]
 [-10.    -9.92  -9.84  -9.76  -9.68]]
[[0.   0.   0.   0.   0.  ]
 [0.08 0.08 0.08 0.08 0.08]
 [0.16 0.16 0.16 0.16 0.16]
 [0.24 0.24 0.24 0.24 0.24]
 [0.32 0.32 0.32 0.32 0.32]]
[[-10.     0.  ]
 [ -9.92   0.  ]
 [ -9.84   0.  ]
 [ -9.76   0.  ]
 [ -9.68   0.  ]]

　 $x_1, x_2$ の値の全ての組み合わせ(格子状の点)を np.meshgrid() で作成して、複製された軸ごとの値(作図用の配列)を x_0_grid, x_1_grid、2つの軸の値をまとめた座標(計算用の配列)を x_arr とします。

　続いて、設定した値に従う確率密度を計算して、グラフを作成していきます。

確率分布の計算

　 $\mathbf{x}$ の値ごとにガウス分布に従う確率密度 $\mathcal{N}(\mathbf{x} \mid \boldsymbol{\mu}, \boldsymbol{\Sigma})$ を計算します。

# 確率密度を計算
dens_grid = multivariate_normal.pdf(
    x=x_arr, mean=mu_d, cov=sigma2_dd
).reshape(x_dims)
print(dens_grid[150:155, 150:155].round(3))

[[0.017 0.017 0.017 0.017 0.016]
 [0.016 0.016 0.016 0.016 0.016]
 [0.016 0.016 0.016 0.015 0.015]
 [0.015 0.015 0.015 0.015 0.015]
 [0.015 0.014 0.014 0.014 0.014]]

　多次元ガウス分布の確率密度は、sp.stats.multivariate_normal.pdf() で計算できます。確率変数の引数 x に $\mathbf{x}$ (の配列)、平均ベクトルの引数 mean に $\boldsymbol{\mu}$ 、分散共分散行列の引数 cov に $\boldsymbol{\Sigma}$ を指定します。

確率分布の作図

　パラメータラベルを作成します。装飾用の処理です。

# ラベル用の文字列を作成
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'$\\mu = ({tmp_mu_str}), '
param_lbl += f'\\Sigma = \\binom{{{tmp_sgm_0_str}}}{{{tmp_sgm_1_str}}}$'
print(param_lbl)

[tex: \mu = (0.0, 10.0), \Sigma = \binom{9.0, 2.5}{2.5, 4.0}$

　ギリシャ文字などの記号を使った数式を表示する場合は、LaTeXの記法(LaTeXコマンド)を使います。数式に変換する範囲(コマンドの文字列)を $ で挟んで、'$コマンド$' の形式で文字列を指定します。ただし、メタ文字(改行文字 \n など)として認識されないために、バックスラッシュを追加(2つ \\ に)する必要がある場合があります。
　オブジェクト(プログラム上の変数)の値を使う場合は、f文字列(や format() )を使います。変数名を波括弧 {} で挟んで、f'{変数名}' の形式で文字列を指定します。(丸め込みなど)数値の文字列を調整する場合は、format() の記法で指定します。
　この例では、行列を表現するのに二項係数 \binom{上段}{下段} で代用しています。
　f文字列の影響を受けずにコマンドの文字列として波括弧を使う場合は、波括弧を二重にする必要があります。この例では、f文字列としても使うため、波括弧を三重にしています。

　ガウス分布の等高線図を作成します。

# ガウス分布を作図
plt.figure(figsize=(9, 6), dpi=100, facecolor='white')
plt.contour(
    x_0_grid, x_1_grid, dens_grid, 
    linewidths=1.0
) # 等高線
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() # カラーバー
plt.grid() # グリッド線
plt.show()

　等高線図は、plt.contour() で描画できます。第1・第2引数に変数、第3引数に確率密度を指定します。

　plt.contourf() を使うと、塗りつぶしの等高線図を作成できます。

# ガウス分布を作図
plt.figure(figsize=(9, 6), dpi=100, facecolor='white')
plt.contourf(
    x_0_grid, x_1_grid, dens_grid
) # 塗りつぶし等高線
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() # カラーバー
plt.grid() # グリッド線
plt.show()

　曲面図を作成します。

# ガウス分布を作図
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=dens_grid, offset=0.0, 
    linewidths=1.0, linestyles=':'
) # 等高線
ax.plot_surface(
    X=x_0_grid, Y=x_1_grid, Z=dens_grid, 
    cmap='viridis', alpha=0.8
) # 曲面
ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$')
ax.set_zlabel('density')
ax.set_title(param_lbl, loc='left')
plt.show()