はじめに
機械学習で登場する確率分布について色々な角度から理解したいシリーズです。
この記事では、Pythonで多次元(多変量)スチューデントのt分布のグラフを作成します。
【前の内容】
【他の記事一覧】
【この記事の内容】
2次元スチューデントのt分布の作図
2次元スチューデントのt分布(Bivariate Student's t Distribution)のグラフを作成します。t分布については「多次元スチューデントのt分布の定義式 - からっぽのしょこ」を参照してください。
利用するパッケージを読み込みます。
# ライブラリを読み込み import numpy as np from scipy.stats import multivariate_t import matplotlib.pyplot as plt from matplotlib.animation import FuncAnimation
定義式の確認
まずは、多次元t分布の定義式を確認します。
t分布は、次の式で定義されます。
$$
\mathrm{St}(\mathbf{x} | \nu, \boldsymbol{\mu}, \boldsymbol{\Sigma})
= \frac{
\Gamma(\frac{\nu + D}{2})
}{
\Gamma(\frac{\nu}{2})
}
\frac{
1
}{
(\pi \nu)^{\frac{D}{2}}
|\boldsymbol{\Sigma}|^{\frac{1}{2}}
}
\left\{
1
+ \frac{1}{\nu}
(\mathbf{x} - \boldsymbol{\mu})^{\top}
\boldsymbol{\Sigma}^{-1}
(\mathbf{x} - \boldsymbol{\mu})
\right\}^{-\frac{\nu+D}{2}}
$$
ここで、$\nu$は形状(自由度)パラメータ、$\boldsymbol{\mu}$は位置ベクトルパラメータ、$\boldsymbol{\Sigma}$はスケール行列パラメータ、$\boldsymbol{\Lambda}$は逆スケール行列パラメータです。また、$\mathbf{A}^{\top}$は行列$\mathbf{A}$の転置行列、$\mathbf{A}^{-1}$は逆行列、$|\mathbf{A}|$は行列式、2分の1乗は平方根$\sqrt{a} = a^{\frac{1}{2}}$です。
$\mathbf{x}, \boldsymbol{\mu}$は$D$次元ベクトル、$\boldsymbol{\Sigma}$は$D \times D$の正定値行列です。
$$
\mathbf{x}
= \begin{pmatrix}
x_1 \\
x_2 \\
\vdots \\
x_D
\end{pmatrix}
,\
\boldsymbol{\mu}
= \begin{pmatrix}
\mu_1 \\
\mu_2 \\
\vdots \\
\mu_D
\end{pmatrix}
,\
\boldsymbol{\Sigma}
= \begin{pmatrix}
\sigma_{1,1} & \sigma_{1,2} & \cdots & \sigma_{1,D} \\
\sigma_{2,1} & \sigma_{2,2} & \cdots & \sigma_{2,D} \\
\vdots & \vdots & \ddots & \vdots \\
\sigma_{D,1} & \sigma_{D,2} & \cdots & \sigma_{D,D}
\end{pmatrix}
$$
$x_d$は実数をとり、$\nu$は正の実数、$\mu_d$は実数、$\sigma_{d,d}$は正の実数、$\sigma_{i,j} = \sigma_{j,i}\ (i \neq j)$は実数、また$\boldsymbol{\Sigma}$は正定値行列を満たす必要があります。
この計算を行いグラフを作成します。
グラフの作成
Matplotlibライブラリを利用して、2次元t分布のグラフを作成します。t分布の確率密度の計算については「【Python】多次元スチューデントのt分布の計算 - からっぽのしょこ」を参照してください。
t分布のパラメータ$\nu, \boldsymbol{\mu}, \boldsymbol{\Sigma}$を設定します。この例では、2次元のグラフで描画するため、次元数を$D = 2$とします。
# 次元数を設定:(固定) D = 2 # 自由度を指定 nu = 3 # 位置ベクトルを指定 mu_d = np.array([6.0, 10.0]) # スケール行列を指定 sigma_dd = np.array( [[1.0, 0.6], [0.6, 4.0]] )
位置ベクトル$\boldsymbol{\mu} = (\mu_1, \mu_2)$、スケール行列$\boldsymbol{\Sigma} = (\sigma_{1,1}, \sigma_{2,1}, \sigma_{1,2}, \sigma_{2,2})$を指定します。
設定したパラメータに応じて、t分布の確率変数がとり得る値$\mathbf{x}$の成分$x_1, x_2$の値を作成します。
# xの値を作成 x_1_vals = np.linspace( start=mu_d[0] - np.sqrt(sigma_dd[0, 0])*3.0, stop=mu_d[0] + np.sqrt(sigma_dd[0, 0])*3.0, num=101 ) x_2_vals = np.linspace( start=mu_d[1] - np.sqrt(sigma_dd[1, 1])*3.0, stop=mu_d[1] + np.sqrt(sigma_dd[1, 1])*3.0, num=101 ) print(x_1_vals[:5]) print(x_2_vals[:5])
[3. 3.06 3.12 3.18 3.24]
[4. 4.12 4.24 4.36 4.48]
$x_1$(x軸)の値をx_1_vals、$x_2$(y軸)の値をx_2_valsとします。この例では、それぞれ$\mu_d$を中心に$\sqrt{\sigma_{d,d}}$の3倍を範囲とします。
$\mathbf{x}$の点を作成します。
# 作図用のxの点を作成 x_1_grid, x_2_grid = np.meshgrid(x_1_vals, x_2_vals) # 計算用のxの点を作成 x_points = np.stack([x_1_grid.flatten(), x_2_grid.flatten()], axis=1) x_dims = x_1_grid.shape print(x_points) print(x_dims)
[[ 3. 4. ]
[ 3.06 4. ]
[ 3.12 4. ]
...
[ 8.88 16. ]
[ 8.94 16. ]
[ 9. 16. ]]
(101, 101)
x_1_valsとx_2_valsの要素の全ての組み合わせ(格子状の点)をnp.meshgrid()で作成してx_1_grid, x_2_gridとします。
x_1_grid, x_2_gridをnp.stack()で列方向に並べてx_pointsとします。x_1_grid, x_2_gridの同じインデックス、x_pointsの各行が点$\mathbf{x} = (x_1, x_2)$に対応します。
$\mathbf{x}$の点ごとの確率密度を計算します。
# 多次元t分布の確率密度を計算 dens = multivariate_t.pdf(x=x_points, loc=mu_d, shape=sigma_dd, df=nu) print(dens)
[0.0011164 0.00116294 0.00121077 ... 0.00121077 0.00116294 0.0011164 ]
多次元t分布の確率密度は、SciPyライブラリのstatsモジュールのmultivariate_tクラスのpdf()で計算できます。確率変数の引数xにx_points、位置ベクトルの引数locにmu_d、スケール行列の引数shapeにsigma_dd、自由度の引数dfにnuを指定します。
パラメータの値を数式で表示するための文字列を作成します。
# パラメータラベル用の文字列を作成 param_text = '$\\nu=' + str(nu) param_text += ', \mu=(' + ', '.join([str(mu) for mu in mu_d]) + ')' param_text += ', \Sigma=' + str([list(sigma_d) for sigma_d in sigma_dd]) + '$' print(param_text)
$\nu=3, \mu=(6.0, 10.0), \Sigma=[[1.0, 0.6], [0.6, 4.0]]$
LaTeXコマンドを使って$で挟みます。
t分布の等高線図を作成します。
# 2次元t分布を作図:(等高線図) plt.figure(figsize=(12, 9)) # 図の設定 cnf = plt.contour(x_1_grid, x_2_grid, dens.reshape(x_dims)) # 等高線 #cnf = plt.contourf(x_1_grid, x_2_grid, dens.reshape(x_dims)) # 塗りつぶし等高線 plt.xlabel('$x_1$') # x軸ラベル plt.ylabel('$x_2$') # y軸ラベル plt.suptitle("Bivariate Student's t-Distribution", fontsize=20) # 全体のタイトル plt.title(param_text, loc='left') # タイトル plt.colorbar(cnf, label='density') # カラーバー plt.grid() # グリッド線 plt.show() # 描画
等高線はplt.contour()、塗りつぶし等高線はplt.contourf()で描画できます。
確率密度の変化を細かく見るには、ヒートマップで可視化します。
# 2次元t分布を作図:(ヒートマップ) plt.figure(figsize=(12, 9)) # 図の設定 pcl = plt.pcolor(x_1_grid, x_2_grid, dens.reshape(x_dims)) # ヒートマップ plt.xlabel('$x_1$') # x軸ラベル plt.ylabel('$x_2$') # y軸ラベル plt.suptitle("Bivariate Student's t-Distribution", fontsize=20) # 全体のタイトル plt.title(param_text, loc='left') # タイトル plt.colorbar(pcl, label='density') # カラーバー plt.grid() # グリッド線 plt.show() # 描画
ヒートマップはplt.pcolor()で描画できます。
3Dプロットを作成します。
# 2次元t分布を作図:曲面図 fig = plt.figure(figsize=(9, 9)) # 図の設定 ax = fig.add_subplot(projection='3d') # 3D用の設定 ax.plot_surface(x_1_grid, x_2_grid, dens.reshape(x_dims), cmap='jet', alpha=0.7) # 曲面図 ax.contour(x_1_grid, x_2_grid, dens.reshape(x_dims), cmap='jet', offset=0.0) # 等高線図 ax.set_xlabel('$x_1$') # x軸ラベル ax.set_ylabel('$x_2$') # y軸ラベル ax.set_zlabel('density') # z軸ラベル fig.suptitle("Bivariate Student's t-Distribution", fontsize=20) # 全体のタイトル ax.set_title(param_text, loc='left') # タイトル plt.show() # 描画
曲面図はax.plot_surface()で描画できます。
t分布のグラフを描画できました。以降は、ここまでの作図処理を用いて、パラメータの影響を確認していきます。
パラメータと分布の形状の関係
パラメータの値を少しずつ変化させて、分布の形状の変化をアニメーションで確認します。
自由度の影響
まずは、自由度$\nu$の値を変化させ、$\boldsymbol{\mu}, \boldsymbol{\Sigma}$を固定します。
・コード(クリックで展開)
# 次元数を設定:(固定) D = 2 # 自由度として利用する値を指定 nu_vals = np.arange(start=0.1, stop=15.1, step=0.1).round(decimals=1) # 位置ベクトルを指定 mu_d = np.array([6.0, 10.0]) # スケール行列を指定 sigma_dd = np.array([[1.0, 0.6], [0.6, 4.0]]) # フレーム数を設定 frame_num = len(nu_vals) print(frame_num)
150
値の間隔が一定になるように$\nu$の値をnu_valsとして作成します。パラメータごとにフレームを切り替えるので、nu_valsの要素数がアニメーションのフレーム数になります。
設定したパラメータに応じて、作図用と計算用の$\mathbf{x}$を作成します。
# xの値を作成 x_1_vals = np.linspace( start=mu_d[0] - np.sqrt(sigma_dd[0, 0])*3.0, stop=mu_d[0] + np.sqrt(sigma_dd[0, 0])*3.0, num=101 ) x_2_vals = np.linspace( start=mu_d[1] - np.sqrt(sigma_dd[1, 1])*3.0, stop=mu_d[1] + np.sqrt(sigma_dd[1, 1])*3.0, num=101 ) # 作図用のxの点を作成 x_1_grid, x_2_grid = np.meshgrid(x_1_vals, x_2_vals) # 計算用のxの点を作成 x_points = np.stack([x_1_grid.flatten(), x_2_grid.flatten()], axis=1) x_shape = x_1_grid.shape print(x_points.round(3))
[[ 3. 4. ]
[ 3.06 4. ]
[ 3.12 4. ]
...
[ 8.88 16. ]
[ 8.94 16. ]
[ 9. 16. ]]
全てのフレームで共通のグラデーションと等高線を引くための値を設定します。
# z軸(確率密度)の最小値・最大値を設定 z_min = 0.0 z_max = multivariate_t.pdf(x=x_points, loc=mu_d, shape=sigma_dd, df=nu_vals.mean()).max() z_max = np.ceil(z_max * 10.0) / 10.0 # 等高線を引く値を設定 z_levels = np.linspace(start=z_min, stop=z_max, num=11) print(z_levels)
[0. 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 ]
z軸(グラデーション)の最小値と最大値を設定して、等高線を描画する値を作成します(もう少し上手く最大値を設定したい。具体的には、小数点以下の桁数を指定して切り上げたい)。
等高線の数は、num引数に対応します。
等高線図のアニメーション(gif画像)を作成します。
# 図を初期化 fig = plt.figure(figsize=(12, 9), facecolor='white') # 図の設定 fig.suptitle("Bivariate Student's t-Distribution", fontsize=20) # 全体のタイトル # カラーバーを表示 tmp = plt.contourf(x_1_grid, x_2_grid, np.zeros(x_shape), vmin=z_min, vmax=z_max, levels=z_levels) # カラーバー用のダミー fig.colorbar(tmp, label='density') # 作図処理を関数として定義 def update(i): # 前フレームのグラフを初期化 plt.cla() # i番目の自由度を取得 nu = nu_vals[i] # 多次元t分布の確率密度を計算 dens = multivariate_t.pdf(x=x_points, loc=mu_d, shape=sigma_dd, df=nu) # パラメータラベル用の文字列を作成 param_text = '$\\nu=' + str(nu) param_text += ', \mu=(' + ', '.join([str(mu) for mu in mu_d]) + ')' param_text += ', \Sigma=' + str([list(sigma_d) for sigma_d in sigma_dd]) + '$' # 2次元t分布を作図:(等高線図) #plt.contour(x_1_grid, x_2_grid, dens.reshape(x_dims), # vmin=z_min, vmax=z_max, levels=z_levels) # 等高線 plt.contourf(x_1_grid, x_2_grid, dens.reshape(x_dims), vmin=z_min, vmax=z_max, levels=z_levels) # 塗りつぶし等高線 plt.xlabel('$x_1$') # x軸ラベル plt.ylabel('$x_2$') # y軸ラベル plt.title(param_text, loc='left') # タイトル plt.grid() # グリッド線 # gif画像を作成 anime_dens = FuncAnimation(fig, update, frames=frame_num, interval=100) # gif画像を保存 anime_dens.save('Bivariate_t_dens_nu_cnf.gif')
各フレームの確率密度の計算と作図処理を関数として定義して、FuncAnimation()を使ってアニメーション(gif画像)を作成する。
同様に、ヒートマップのアニメーションを作成します。
# 図を初期化 fig = plt.figure(figsize=(12, 9), facecolor='white') # 図の設定 fig.suptitle("Bivariate Student's t-Distribution", fontsize=20) # 全体のタイトル # カラーバーを表示 tmp = plt.pcolor(x_1_grid, x_2_grid, np.zeros(x_shape), vmin=z_min, vmax=z_max) # カラーバー用のダミー fig.colorbar(tmp, label='density') # 作図処理を関数として定義 def update(i): # 前フレームのグラフを初期化 plt.cla() # i番目の自由度を取得 nu = nu_vals[i] # 多次元t分布の確率密度を計算 dens = multivariate_t.pdf(x=x_points, loc=mu_d, shape=sigma_dd, df=nu) # パラメータラベル用の文字列を作成 param_text = '$\\nu=' + str(nu) param_text += ', \mu=(' + ', '.join([str(mu) for mu in mu_d]) + ')' param_text += ', \Sigma=' + str([list(sigma_d) for sigma_d in sigma_dd]) + '$' # 2次元t分布を作図:(ヒートマップ) plt.pcolor(x_1_grid, x_2_grid, dens.reshape(x_dims), vmin=z_min, vmax=z_max) # 塗りつぶし等高線 plt.xlabel('$x_1$') # x軸ラベル plt.ylabel('$x_2$') # y軸ラベル plt.title(param_text, loc='left') # タイトル plt.grid() # グリッド線 # gif画像を作成 anime_dens = FuncAnimation(fig, update, frames=frame_num, interval=100) # gif画像を保存 anime_dens.save('Bivariate_t_dens_nu_pcl.gif')
3Dプロットのアニメーションを作成します。
# 図を初期化 fig = plt.figure(figsize=(9, 9), facecolor='white') # 図の設定 ax = fig.add_subplot(projection='3d') # 3D用の設定 fig.suptitle("Bivariate Student's t-Distribution", fontsize=20) # 全体のタイトル # 作図処理を関数として定義 def update(i): # 前フレームのグラフを初期化 plt.cla() # i番目の自由度を取得 nu = nu_vals[i] # 多次元t分布の確率密度を計算 dens = multivariate_t.pdf(x=x_points, loc=mu_d, shape=sigma_dd, df=nu) # パラメータラベル用の文字列を作成 param_text = '$\\nu=' + str(nu) param_text += ', \mu=(' + ', '.join([str(mu) for mu in mu_d]) + ')' param_text += ', \Sigma=' + str([list(sigma_d) for sigma_d in sigma_dd]) + '$' # 2次元t分布を作図:(曲面図) ax.plot_surface(x_1_grid, x_2_grid, dens.reshape(x_shape), cmap='jet', vmin=z_min, vmax=z_max, alpha=0.8) # 曲面図 ax.contour(x_1_grid, x_2_grid, dens.reshape(x_shape), cmap='jet', vmin=z_min, vmax=z_max, levels=z_levels, offset=0.0) # 等高線図 ax.set_xlabel('$x_1$') # x軸ラベル ax.set_ylabel('$x_2$') # y軸ラベル ax.set_zlabel('density') # z軸ラベル ax.set_title(param_text, loc='left') # タイトル # gif画像を作成 anime_dens = FuncAnimation(fig, update, frames=frame_num, interval=100) # gif画像を保存 anime_dens.save('Bivariate_t_dens_nu_srf.gif')
$\nu$の値に応じて、形状が変化するのが分かります。
位置ベクトル(1軸)の影響
次は、x軸方向の位置パラメータ$\mu_1$の値を変化させ、$\nu, \mu_2, \boldsymbol{\Sigma}$を固定します。
・コード(クリックで展開)
# 次元数を設定:(固定) D = 2 # 自由度を指定 nu = 3 # x軸の位置パラメータとして利用する値を指定 mu_1_vals = np.linspace(start=-2.0, stop=2.0, num=101).round(decimals=2) # y軸位置パラメータを指定 mu_2 = 10.0 # スケール行列を指定 sigma_dd = np.array([[1.0, 0.6], [0.6, 4.0]]) # フレーム数を設定 frame_num = len(mu_1_vals) print(frame_num)
101
値の間隔が一定になるように$\mu_1$の値をmu_1_valsとして作成します。また、$\mu_2$をmu_2として値を指定します。
作図用と計算用の$\mathbf{x}$を作成します。
# xの値を作成 x_1_vals = np.linspace( start=mu_1_vals.min() - np.sqrt(sigma_dd[0, 0])*2.0, stop=mu_1_vals.max() + np.sqrt(sigma_dd[0, 0])*2.0, num=101 ) x_2_vals = np.linspace( start=mu_2 - np.sqrt(sigma_dd[1, 1])*3.0, stop=mu_2 + np.sqrt(sigma_dd[1, 1])*3.0, num=101 ) # 作図用のxの点を作成 x_1_grid, x_2_grid = np.meshgrid(x_1_vals, x_2_vals) # 計算用のxの点を作成 x_points = np.stack([x_1_grid.flatten(), x_2_grid.flatten()], axis=1) x_shape = x_1_grid.shape print(x_points.round(3))
[[-4. 4. ]
[-3.92 4. ]
[-3.84 4. ]
...
[ 3.84 16. ]
[ 3.92 16. ]
[ 4. 16. ]]
等高線の設定用の配列を作成します。
# z軸(確率密度)の最小値・最大値を設定 z_min = 0.0 z_max = multivariate_t.pdf(x=x_points, loc=[mu_1_vals.mean(), mu_2], shape=sigma_dd, df=nu).max() z_max = np.ceil(z_max * 10.0) / 10.0 # 等高線を引く値を設定 z_levels = np.linspace(start=z_min, stop=z_max, num=11) print(z_levels)
[0. 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 ]
等高線図のアニメーションを作成します。
# 図を初期化 fig = plt.figure(figsize=(12, 9), facecolor='white') # 図の設定 fig.suptitle("Bivariate Student's t-Distribution", fontsize=20) # 全体のタイトル # カラーバーを表示 tmp = plt.contourf(x_1_grid, x_2_grid, np.zeros(x_shape), vmin=z_min, vmax=z_max, levels=z_levels) # カラーバー用のダミー fig.colorbar(tmp, label='density') # 作図処理を関数として定義 def update(i): # 前フレームのグラフを初期化 plt.cla() # i番目の位置ベクトルを取得 mu_1 = mu_1_vals[i] mu_d = np.array([mu_1, mu_2]) # 多次元t分布の確率密度を計算 dens = multivariate_t.pdf(x=x_points, loc=mu_d, shape=sigma_dd, df=nu) # パラメータラベル用の文字列を作成 param_text = '$\\nu=' + str(nu) param_text += ', \mu=(' + ', '.join([str(mu) for mu in mu_d]) + ')' param_text += ', \Sigma=' + str([list(sigma_d) for sigma_d in sigma_dd]) + '$' # 2次元t分布を作図:(等高線図) #plt.contour(x_1_grid, x_2_grid, dens.reshape(x_dims), # vmin=z_min, vmax=z_max, levels=z_levels) # 等高線 plt.contourf(x_1_grid, x_2_grid, dens.reshape(x_dims), vmin=z_min, vmax=z_max, levels=z_levels) # 塗りつぶし等高線 plt.xlabel('$x_1$') # x軸ラベル plt.ylabel('$x_2$') # y軸ラベル plt.title(param_text, loc='left') # タイトル plt.grid() # グリッド線 # gif画像を作成 anime_dens = FuncAnimation(fig, update, frames=frame_num, interval=100) # gif画像を保存 anime_dens.save('Bivariate_t_dens_mu1_cnf.gif')
ヒートマップのアニメーションを作成します。
# 図を初期化 fig = plt.figure(figsize=(12, 9), facecolor='white') # 図の設定 fig.suptitle("Bivariate Student's t-Distribution", fontsize=20) # 全体のタイトル # カラーバーを表示 tmp = plt.pcolor(x_1_grid, x_2_grid, np.zeros(x_shape), vmin=z_min, vmax=z_max) # カラーバー用のダミー fig.colorbar(tmp, label='density') # 作図処理を関数として定義 def update(i): # 前フレームのグラフを初期化 plt.cla() # i番目の位置ベクトルを取得 mu_1 = mu_1_vals[i] mu_d = np.array([mu_1, mu_2]) # 多次元t分布の確率密度を計算 dens = multivariate_t.pdf(x=x_points, loc=mu_d, shape=sigma_dd, df=nu) # パラメータラベル用の文字列を作成 param_text = '$\\nu=' + str(nu) param_text += ', \mu=(' + ', '.join([str(mu) for mu in mu_d]) + ')' param_text += ', \Sigma=' + str([list(sigma_d) for sigma_d in sigma_dd]) + '$' # 2次元t分布を作図:(ヒートマップ) plt.pcolor(x_1_grid, x_2_grid, dens.reshape(x_dims), vmin=z_min, vmax=z_max) # 塗りつぶし等高線 plt.xlabel('$x_1$') # x軸ラベル plt.ylabel('$x_2$') # y軸ラベル plt.title(param_text, loc='left') # タイトル plt.grid() # グリッド線 # gif画像を作成 anime_dens = FuncAnimation(fig, update, frames=frame_num, interval=100) # gif画像を保存 anime_dens.save('Bivariate_t_dens_mu1_pcl.gif')
曲面図のアニメーションを作成します。
# 図を初期化 fig = plt.figure(figsize=(9, 9), facecolor='white') # 図の設定 ax = fig.add_subplot(projection='3d') # 3D用の設定 fig.suptitle("Bivariate Student's t-Distribution", fontsize=20) # 全体のタイトル # 作図処理を関数として定義 def update(i): # 前フレームのグラフを初期化 plt.cla() # i番目のパラメータを取得 mu_1 = mu_1_vals[i] mu_d = np.array([mu_1, mu_2]) # 多次元t分布の確率密度を計算 dens = multivariate_t.pdf(x=x_points, loc=mu_d, shape=sigma_dd, df=nu) # パラメータラベル用の文字列を作成 param_text = '$\\nu=' + str(nu) param_text += ', \mu=(' + ', '.join([str(mu) for mu in mu_d]) + ')' param_text += ', \Sigma=' + str([list(sigma_d) for sigma_d in sigma_dd]) + '$' # 2次元t分布を作図:(曲面図) ax.plot_surface(x_1_grid, x_2_grid, dens.reshape(x_dims), cmap='jet', vmin=z_min, vmax=z_max, alpha=0.8) # 曲面図 ax.contour(x_1_grid, x_2_grid, dens.reshape(x_dims), cmap='jet', vmin=z_min, vmax=z_max, levels=z_levels, offset=0.0) # 等高線図 ax.set_xlabel('$x_1$') # x軸ラベル ax.set_ylabel('$x_2$') # y軸ラベル ax.set_zlabel('density') # z軸ラベル ax.set_title(param_text, loc='left') # タイトル ax.set_zlim(zmin=z_min, zmax=z_max) # z軸の表示範囲 # gif画像を作成 anime_dens = FuncAnimation(fig, update, frames=frame_num, interval=100) # gif画像を保存 anime_dens.save('Bivariate_t_dens_mu1_srf.gif')
$\mu_1$の値に応じて、x軸方向に移動するのが分かります。
位置ベクトル(2軸)の影響
y軸方向の位置パラメータ$\mu_2$の値を変化させ、$\nu, \mu_1, \boldsymbol{\Sigma}$を固定します。
・コード(クリックで展開)
# 次元数を設定:(固定) D = 2 # 自由度を指定 nu = 3 # x軸位置パラメータを指定 mu_1 = 6.0 # y軸の位置パラメータとして利用する値を指定 mu_2_vals = np.linspace(start=-2.0, stop=2.0, num=101).round(decimals=2) # スケール行列を指定 sigma_dd = np.array([[1.0, 0.6], [0.6, 4.0]]) # フレーム数を設定 frame_num = len(mu_2_vals) print(frame_num)
101
値の間隔が一定になるように$\mu_2$の値をmu_2_valsとして作成します。また、$\mu_1$をmu_1として値を指定します。
作図用と計算用の$\mathbf{x}$を作成します。
# xの値を作成 x_1_vals = np.linspace( start=mu_1 - np.sqrt(sigma_dd[0, 0])*3.0, stop=mu_1 + np.sqrt(sigma_dd[0, 0])*3.0, num=101 ) x_2_vals = np.linspace( start=mu_2_vals.min() - np.sqrt(sigma_dd[1, 1])*2.0, stop=mu_2_vals.max() + np.sqrt(sigma_dd[1, 1])*2.0, num=101 ) # 作図用のxの点を作成 x_1_grid, x_2_grid = np.meshgrid(x_1_vals, x_2_vals) # 計算用のxの点を作成 x_points = np.stack([x_1_grid.flatten(), x_2_grid.flatten()], axis=1) x_shape = x_1_grid.shape print(x_points.round(3))
[[ 3. -6. ]
[ 3.06 -6. ]
[ 3.12 -6. ]
...
[ 8.88 6. ]
[ 8.94 6. ]
[ 9. 6. ]]
等高線の設定用の配列を作成します。
# z軸(確率密度)の最小値・最大値を設定 z_min = 0.0 z_max = multivariate_t.pdf(x=x_points, loc=[mu_1, mu_2_vals.mean()], shape=sigma_dd, df=nu).max() z_max = np.ceil(z_max * 10.0) / 10.0 # 等高線を引く値を設定 z_levels = np.linspace(start=z_min, stop=z_max, num=11) print(z_levels)
[0. 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 ]
等高線図のアニメーションを作成します。
# 図を初期化 fig = plt.figure(figsize=(12, 9), facecolor='white') # 図の設定 fig.suptitle("Bivariate Student's t-Distribution", fontsize=20) # 全体のタイトル # カラーバーを表示 tmp = plt.contourf(x_1_grid, x_2_grid, np.zeros(x_shape), vmin=z_min, vmax=z_max, levels=z_levels) # カラーバー用のダミー fig.colorbar(tmp, label='density') # 作図処理を関数として定義 def update(i): # 前フレームのグラフを初期化 plt.cla() # i番目の位置ベクトルを取得 mu_2 = mu_2_vals[i] mu_d = np.array([mu_1, mu_2]) # 多次元t分布の確率密度を計算 dens = multivariate_t.pdf(x=x_points, loc=mu_d, shape=sigma_dd, df=nu) # パラメータラベル用の文字列を作成 param_text = '$\\nu=' + str(nu) param_text += ', \mu=(' + ', '.join([str(mu) for mu in mu_d]) + ')' param_text += ', \Sigma=' + str([list(sigma_d) for sigma_d in sigma_dd]) + '$' # 2次元t分布を作図:(等高線図) #plt.contour(x_1_grid, x_2_grid, dens.reshape(x_dims), # vmin=z_min, vmax=z_max, levels=z_levels) # 等高線 plt.contourf(x_1_grid, x_2_grid, dens.reshape(x_dims), vmin=z_min, vmax=z_max, levels=z_levels) # 塗りつぶし等高線 plt.xlabel('$x_1$') # x軸ラベル plt.ylabel('$x_2$') # y軸ラベル plt.title(param_text, loc='left') # タイトル plt.grid() # グリッド線 # gif画像を作成 anime_dens = FuncAnimation(fig, update, frames=frame_num, interval=100) # gif画像を保存 anime_dens.save('Bivariate_t_dens_mu2_cnf.gif')
ヒートマップのアニメーションを作成します。
# 図を初期化 fig = plt.figure(figsize=(12, 9), facecolor='white') # 図の設定 fig.suptitle("Bivariate Student's t-Distribution", fontsize=20) # 全体のタイトル # カラーバーを表示 tmp = plt.pcolor(x_1_grid, x_2_grid, np.zeros(x_shape), vmin=z_min, vmax=z_max) # カラーバー用のダミー fig.colorbar(tmp, label='density') # 作図処理を関数として定義 def update(i): # 前フレームのグラフを初期化 plt.cla() # i番目の位置ベクトルを取得 mu_2 = mu_2_vals[i] mu_d = np.array([mu_1, mu_2]) # 多次元t分布の確率密度を計算 dens = multivariate_t.pdf(x=x_points, loc=mu_d, shape=sigma_dd, df=nu) # パラメータラベル用の文字列を作成 param_text = '$\\nu=' + str(nu) param_text += ', \mu=(' + ', '.join([str(mu) for mu in mu_d]) + ')' param_text += ', \Sigma=' + str([list(sigma_d) for sigma_d in sigma_dd]) + '$' # 2次元t分布を作図:(ヒートマップ) plt.pcolor(x_1_grid, x_2_grid, dens.reshape(x_dims), vmin=z_min, vmax=z_max) # 塗りつぶし等高線 plt.xlabel('$x_1$') # x軸ラベル plt.ylabel('$x_2$') # y軸ラベル plt.title(param_text, loc='left') # タイトル plt.grid() # グリッド線 # gif画像を作成 anime_dens = FuncAnimation(fig, update, frames=frame_num, interval=100) # gif画像を保存 anime_dens.save('Bivariate_t_dens_mu2_pcl.gif')
曲面図のアニメーションを作成します。
# 図を初期化 fig = plt.figure(figsize=(9, 9), facecolor='white') # 図の設定 ax = fig.add_subplot(projection='3d') # 3D用の設定 fig.suptitle("Bivariate Student's t-Distribution", fontsize=20) # 全体のタイトル # 作図処理を関数として定義 def update(i): # 前フレームのグラフを初期化 plt.cla() # i番目のパラメータを取得 mu_2 = mu_2_vals[i] mu_d = np.array([mu_1, mu_2]) # 多次元t分布の確率密度を計算 dens = multivariate_t.pdf(x=x_points, loc=mu_d, shape=sigma_dd, df=nu) # パラメータラベル用の文字列を作成 param_text = '$\\nu=' + str(nu) param_text += ', \mu=(' + ', '.join([str(mu) for mu in mu_d]) + ')' param_text += ', \Sigma=' + str([list(sigma_d) for sigma_d in sigma_dd]) + '$' # 2次元t分布を作図:(曲面図) ax.plot_surface(x_1_grid, x_2_grid, dens.reshape(x_dims), cmap='jet', vmin=z_min, vmax=z_max, alpha=0.8) # 曲面図 ax.contour(x_1_grid, x_2_grid, dens.reshape(x_dims), cmap='jet', vmin=z_min, vmax=z_max, levels=z_levels, offset=0.0) # 等高線図 ax.set_xlabel('$x_1$') # x軸ラベル ax.set_ylabel('$x_2$') # y軸ラベル ax.set_zlabel('density') # z軸ラベル ax.set_title(param_text, loc='left') # タイトル ax.set_zlim(zmin=z_min, zmax=z_max) # z軸の表示範囲 # gif画像を作成 anime_dens = FuncAnimation(fig, update, frames=frame_num, interval=100) # gif画像を保存 anime_dens.save('Bivariate_t_dens_mu2_srf.gif')
$\mu_2$の値に応じて、y軸方向に移動するのが分かります。
スケール行列(1,1成分)の影響
続いて、x軸方向のスケール$\sigma_{1,1}$の値を変化させ、$\nu, \boldsymbol{\mu}, \sigma_{1,2}, \sigma_{2,1}, \sigma_{2,2}$を固定します。
・コード(クリックで展開)
# 次元数を設定:(固定) D = 2 # 自由度を指定 nu = 3 # 位置ベクトルを指定 mu_1 = np.array([6.0, 10.0]) # x軸のスケールパラメータとして利用する値を指定 sigma_11_vals = np.arange(start=0.5, stop=6.0, step=0.1).round(decimals=1) # y軸のスケールパラメータを指定 sigma_22 = 4.0 # x・y軸のスケールパラメータを指定 sigma_12 = 0.6 # フレーム数を設定 frame_num = len(sigma_11_vals) print(frame_num)
55
値の間隔が一定になるように$\sigma_{1,1}$の値をsigma_11_valsとして作成します。また、$\sigma_{2,2}$をsigma_22、$\sigma_{1,2} = \sigma_{2,1}$をsigma_12として値を指定します。
作図用と計算用の$\mathbf{x}$を作成します。
# xの値を作成 x_1_vals = np.linspace( start=mu_d[0] - np.sqrt(sigma_11_vals.max())*2.0, stop=mu_d[0] + np.sqrt(sigma_11_vals.max())*2.0, num=101 ) x_2_vals = np.linspace( start=mu_d[1] - np.sqrt(sigma_22)*3.0, stop=mu_d[1] + np.sqrt(sigma_22)*3.0, num=101 ) # 作図用のxの点を作成 x_1_grid, x_2_grid = np.meshgrid(x_1_vals, x_2_vals) # 計算用のxの点を作成 x_points = np.stack([x_1_grid.flatten(), x_2_grid.flatten()], axis=1) x_shape = x_1_grid.shape print(x_points.round(3))
[[ 1.142 4. ]
[ 1.239 4. ]
[ 1.336 4. ]
...
[10.664 16. ]
[10.761 16. ]
[10.858 16. ]]
等高線の設定用の配列を作成します。
# z軸(確率密度)の最小値・最大値を設定 z_min = 0.0 z_max = multivariate_t.pdf(x=x_points, loc=mu_d, shape=np.array([[sigma_11_vals.min(), sigma_12], [sigma_12, sigma_22]]), df=nu).max() z_max = np.ceil(z_max * 10.0) / 10.0 # 等高線を引く値を設定 z_levels = np.linspace(start=z_min, stop=z_max, num=11) print(z_levels)
[0. 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 ]
等高線図のアニメーションを作成します。
# 図を初期化 fig = plt.figure(figsize=(12, 9), facecolor='white') # 図の設定 fig.suptitle("Bivariate Student's t-Distribution", fontsize=20) # 全体のタイトル # カラーバーを表示 tmp = plt.contourf(x_1_grid, x_2_grid, np.zeros(x_shape), vmin=z_min, vmax=z_max, levels=z_levels) # カラーバー用のダミー fig.colorbar(tmp, label='density') # 作図処理を関数として定義 def update(i): # 前フレームのグラフを初期化 plt.cla() # i番目のスケール行列を取得 sigma_11 = sigma_11_vals[i] sigma_dd = np.array([[sigma_11, sigma_12], [sigma_12, sigma_22]]) # 多次元t分布の確率密度を計算 dens = multivariate_t.pdf(x=x_points, loc=mu_d, shape=sigma_dd, df=nu) # パラメータラベル用の文字列を作成 param_text = '$\\nu=' + str(nu) param_text += ', \mu=(' + ', '.join([str(mu) for mu in mu_d]) + ')' param_text += ', \Sigma=' + str([list(sigma_d) for sigma_d in sigma_dd]) + '$' # 2次元t分布を作図:(等高線図) #plt.contour(x_1_grid, x_2_grid, dens.reshape(x_dims), # vmin=z_min, vmax=z_max, levels=z_levels) # 等高線 plt.contourf(x_1_grid, x_2_grid, dens.reshape(x_dims), vmin=z_min, vmax=z_max, levels=z_levels) # 塗りつぶし等高線 plt.xlabel('$x_1$') # x軸ラベル plt.ylabel('$x_2$') # y軸ラベル plt.title(param_text, loc='left') # タイトル plt.grid() # グリッド線 # gif画像を作成 anime_dens = FuncAnimation(fig, update, frames=frame_num, interval=100) # gif画像を保存 anime_dens.save('Bivariate_t_dens_sigma11_cnf.gif')
ヒートマップのアニメーションを作成します。
# 図を初期化 fig = plt.figure(figsize=(12, 9), facecolor='white') # 図の設定 fig.suptitle("Bivariate Student's t-Distribution", fontsize=20) # 全体のタイトル # カラーバーを表示 tmp = plt.pcolor(x_1_grid, x_2_grid, np.zeros(x_shape), vmin=z_min, vmax=z_max) # カラーバー用のダミー fig.colorbar(tmp, label='density') # 作図処理を関数として定義 def update(i): # 前フレームのグラフを初期化 plt.cla() # i番目のスケール行列を取得 sigma_11 = sigma_11_vals[i] sigma_dd = np.array([[sigma_11, sigma_12], [sigma_12, sigma_22]]) # 多次元t分布の確率密度を計算 dens = multivariate_t.pdf(x=x_points, loc=mu_d, shape=sigma_dd, df=nu) # パラメータラベル用の文字列を作成 param_text = '$\\nu=' + str(nu) param_text += ', \mu=(' + ', '.join([str(mu) for mu in mu_d]) + ')' param_text += ', \Sigma=' + str([list(sigma_d) for sigma_d in sigma_dd]) + '$' # 2次元t分布を作図:(ヒートマップ) plt.pcolor(x_1_grid, x_2_grid, dens.reshape(x_dims), vmin=z_min, vmax=z_max) # 塗りつぶし等高線 plt.xlabel('$x_1$') # x軸ラベル plt.ylabel('$x_2$') # y軸ラベル plt.title(param_text, loc='left') # タイトル plt.grid() # グリッド線 # gif画像を作成 anime_dens = FuncAnimation(fig, update, frames=frame_num, interval=100) # gif画像を保存 anime_dens.save('Bivariate_t_dens_sigma11_pcl.gif')
曲面図のアニメーションを作成します。
# 図を初期化 fig = plt.figure(figsize=(9, 9), facecolor='white') # 図の設定 ax = fig.add_subplot(projection='3d') # 3D用の設定 fig.suptitle("Bivariate Student's t-Distribution", fontsize=20) # 全体のタイトル # 作図処理を関数として定義 def update(i): # 前フレームのグラフを初期化 plt.cla() # i番目のスケール行列を取得 sigma_11 = sigma_11_vals[i] sigma_dd = np.array([[sigma_11, sigma_12], [sigma_12, sigma_22]]) # 多次元t分布の確率密度を計算 dens = multivariate_t.pdf(x=x_points, loc=mu_d, shape=sigma_dd, df=nu) # パラメータラベル用の文字列を作成 param_text = '$\\nu=' + str(nu) param_text += ', \mu=(' + ', '.join([str(mu) for mu in mu_d]) + ')' param_text += ', \Sigma=' + str([list(sigma_d) for sigma_d in sigma_dd]) + '$' # 2次元t分布を作図:(曲面図) ax.plot_surface(x_1_grid, x_2_grid, dens.reshape(x_dims), cmap='jet', vmin=z_min, vmax=z_max, alpha=0.8) # 曲面図 ax.contour(x_1_grid, x_2_grid, dens.reshape(x_dims), cmap='jet', vmin=z_min, vmax=z_max, levels=z_levels, offset=0.0) # 等高線図 ax.set_xlabel('$x_1$') # x軸ラベル ax.set_ylabel('$x_2$') # y軸ラベル ax.set_zlabel('density') # z軸ラベル ax.set_title(param_text, loc='left') # タイトル ax.set_zlim(zmin=z_min, zmax=z_max) # z軸の表示範囲 # gif画像を作成 anime_dens = FuncAnimation(fig, update, frames=frame_num, interval=100) # gif画像を保存 anime_dens.save('Bivariate_t_dens_sigma11_srf.gif')
$\sigma_{1,1}$の値に応じて、x軸方向に広がるのが分かります。
スケール行列(2,2成分)の影響
y軸方向のスケール$\sigma_{2,2}$の値を変化させ、$\nu, \boldsymbol{\mu}, \sigma_{1,1}, \sigma_{1,2}, \sigma_{2,1}$を固定します。
・コード(クリックで展開)
# 次元数を設定:(固定) D = 2 # 自由度を指定 nu = 3 # 位置ベクトルを指定 mu_1 = np.array([6.0, 10.0]) # x軸のスケールパラメータを指定 sigma_11 = 1.0 # y軸のスケールパラメータとして利用する値を指定 sigma_22_vals = np.arange(start=0.5, stop=6.0, step=0.1).round(decimals=1) # x・y軸のスケールパラメータを指定 sigma_12 = 0.6 # フレーム数を設定 frame_num = len(sigma_22_vals) print(frame_num)
55
値の間隔が一定になるように$\sigma_{2,2}$の値をsigma_22_valsとして作成します。また、$\sigma_{1,1}$をsigma_11、$\sigma_{1,2} = \sigma_{2,1}$をsigma_12として値を指定します。
作図用と計算用の$\mathbf{x}$を作成します。
# xの値を作成 x_1_vals = np.linspace( start=mu_d[0] - np.sqrt(sigma_11)*3.0, stop=mu_d[0] + np.sqrt(sigma_11)*3.0, num=101 ) x_2_vals = np.linspace( start=mu_d[1] - np.sqrt(sigma_22_vals.max())*2.0, stop=mu_d[1] + np.sqrt(sigma_22_vals.max())*2.0, num=101 ) # 作図用のxの点を作成 x_1_grid, x_2_grid = np.meshgrid(x_1_vals, x_2_vals) # 計算用のxの点を作成 x_points = np.stack([x_1_grid.flatten(), x_2_grid.flatten()], axis=1) x_shape = x_1_grid.shape print(x_points.round(3))
[[ 3. 5.142]
[ 3.06 5.142]
[ 3.12 5.142]
...
[ 8.88 14.858]
[ 8.94 14.858]
[ 9. 14.858]]
等高線の設定用の配列を作成します。
# z軸(確率密度)の最小値・最大値を設定 z_min = 0.0 z_max = multivariate_t.pdf(x=x_points, loc=mu_d, shape=np.array([[sigma_11, sigma_12], [sigma_12, sigma_22_vals.min()]]), df=nu).max() z_max = np.ceil(z_max * 10.0) / 10.0 # 等高線を引く値を設定 z_levels = np.linspace(start=z_min, stop=z_max, num=11) print(z_levels)
[0. 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 ]
等高線図のアニメーションを作成します。
# 図を初期化 fig = plt.figure(figsize=(12, 9), facecolor='white') # 図の設定 fig.suptitle("Bivariate Student's t-Distribution", fontsize=20) # 全体のタイトル # カラーバーを表示 tmp = plt.contourf(x_1_grid, x_2_grid, np.zeros(x_shape), vmin=z_min, vmax=z_max, levels=z_levels) # カラーバー用のダミー fig.colorbar(tmp, label='density') # 作図処理を関数として定義 def update(i): # 前フレームのグラフを初期化 plt.cla() # i番目のスケール行列を取得 sigma_22 = sigma_22_vals[i] sigma_dd = np.array([[sigma_11, sigma_12], [sigma_12, sigma_22]]) # 多次元t分布の確率密度を計算 dens = multivariate_t.pdf(x=x_points, loc=mu_d, shape=sigma_dd, df=nu) # パラメータラベル用の文字列を作成 param_text = '$\\nu=' + str(nu) param_text += ', \mu=(' + ', '.join([str(mu) for mu in mu_d]) + ')' param_text += ', \Sigma=' + str([list(sigma_d) for sigma_d in sigma_dd]) + '$' # 2次元t分布を作図:(等高線図) #plt.contour(x_1_grid, x_2_grid, dens.reshape(x_dims), # vmin=z_min, vmax=z_max, levels=z_levels) # 等高線 plt.contourf(x_1_grid, x_2_grid, dens.reshape(x_dims), vmin=z_min, vmax=z_max, levels=z_levels) # 塗りつぶし等高線 plt.xlabel('$x_1$') # x軸ラベル plt.ylabel('$x_2$') # y軸ラベル plt.title(param_text, loc='left') # タイトル plt.grid() # グリッド線 # gif画像を作成 anime_dens = FuncAnimation(fig, update, frames=frame_num, interval=100) # gif画像を保存 anime_dens.save('Bivariate_t_dens_sigma22_cnf.gif')
ヒートマップのアニメーションを作成します。
# 図を初期化 fig = plt.figure(figsize=(12, 9), facecolor='white') # 図の設定 fig.suptitle("Bivariate Student's t-Distribution", fontsize=20) # 全体のタイトル # カラーバーを表示 tmp = plt.pcolor(x_1_grid, x_2_grid, np.zeros(x_shape), vmin=z_min, vmax=z_max) # カラーバー用のダミー fig.colorbar(tmp, label='density') # 作図処理を関数として定義 def update(i): # 前フレームのグラフを初期化 plt.cla() # i番目のスケール行列を取得 sigma_22 = sigma_22_vals[i] sigma_dd = np.array([[sigma_11, sigma_12], [sigma_12, sigma_22]]) # 多次元t分布の確率密度を計算 dens = multivariate_t.pdf(x=x_points, loc=mu_d, shape=sigma_dd, df=nu) # パラメータラベル用の文字列を作成 param_text = '$\\nu=' + str(nu) param_text += ', \mu=(' + ', '.join([str(mu) for mu in mu_d]) + ')' param_text += ', \Sigma=' + str([list(sigma_d) for sigma_d in sigma_dd]) + '$' # 2次元t分布を作図:(ヒートマップ) plt.pcolor(x_1_grid, x_2_grid, dens.reshape(x_dims), vmin=z_min, vmax=z_max) # 塗りつぶし等高線 plt.xlabel('$x_1$') # x軸ラベル plt.ylabel('$x_2$') # y軸ラベル plt.title(param_text, loc='left') # タイトル plt.grid() # グリッド線 # gif画像を作成 anime_dens = FuncAnimation(fig, update, frames=frame_num, interval=100) # gif画像を保存 anime_dens.save('Bivariate_t_dens_sigma22_pcl.gif')
曲面図のアニメーションを作成します。
# 図を初期化 fig = plt.figure(figsize=(9, 9), facecolor='white') # 図の設定 ax = fig.add_subplot(projection='3d') # 3D用の設定 fig.suptitle("Bivariate Student's t-Distribution", fontsize=20) # 全体のタイトル # 作図処理を関数として定義 def update(i): # 前フレームのグラフを初期化 plt.cla() # i番目のスケール行列を取得 sigma_22 = sigma_22_vals[i] sigma_dd = np.array([[sigma_11, sigma_12], [sigma_12, sigma_22]]) # 多次元t分布の確率密度を計算 dens = multivariate_t.pdf(x=x_points, loc=mu_d, shape=sigma_dd, df=nu) # パラメータラベル用の文字列を作成 param_text = '$\\nu=' + str(nu) param_text += ', \mu=(' + ', '.join([str(mu) for mu in mu_d]) + ')' param_text += ', \Sigma=' + str([list(sigma_d) for sigma_d in sigma_dd]) + '$' # 2次元t分布を作図:(曲面図) ax.plot_surface(x_1_grid, x_2_grid, dens.reshape(x_dims), cmap='jet', vmin=z_min, vmax=z_max, alpha=0.8) # 曲面図 ax.contour(x_1_grid, x_2_grid, dens.reshape(x_dims), cmap='jet', vmin=z_min, vmax=z_max, levels=z_levels, offset=0.0) # 等高線図 ax.set_xlabel('$x_1$') # x軸ラベル ax.set_ylabel('$x_2$') # y軸ラベル ax.set_zlabel('density') # z軸ラベル ax.set_title(param_text, loc='left') # タイトル ax.set_zlim(zmin=z_min, zmax=z_max) # z軸の表示範囲 # gif画像を作成 anime_dens = FuncAnimation(fig, update, frames=frame_num, interval=100) # gif画像を保存 anime_dens.save('Bivariate_t_dens_sigma22_srf.gif')
$\sigma_{2,2}$の値に応じて、y軸方向に広がるのが分かります。
スケール行列(1,2成分)の影響
両軸のスケール$\sigma_{1,2}, \sigma_{2,1}$の値を変化させ、$\nu, \boldsymbol{\mu}, \sigma_{1,1}, \sigma_{2,2}$を固定します。
・コード(クリックで展開)
# 次元数を設定:(固定) D = 2 # 自由度を指定 nu = 3 # 位置ベクトルを指定 mu_1 = np.array([6.0, 10.0]) # x軸・y軸のスケールパラメータを指定 sigma_11 = 4.0 sigma_22 = 10.0 # x・y軸のスケールパラメータとして利用する値を指定 sigma_12_vals = np.linspace(start=-5.0, stop=5.0, num=101).round(decimals=1) # フレーム数を設定 frame_num = len(sigma_12_vals) print(frame_num)
101
値の間隔が一定になるように$\sigma_{1,2} = \sigma_{2,1}$の値をsigma_12_valsとして作成します。また、$\sigma_{1,1}, \sigma_{2,2}$をsigma_11, sigma_22として値を指定します。
作図用と計算用の$\mathbf{x}$を作成します。
# xの値を作成 x_1_vals = np.linspace( start=mu_d[0] - np.sqrt(sigma_11)*3.0, stop=mu_d[0] + np.sqrt(sigma_11)*3.0, num=101 ) x_2_vals = np.linspace( start=mu_d[1] - np.sqrt(sigma_22)*3.0, stop=mu_d[1] + np.sqrt(sigma_22)*3.0, num=101 ) # 作図用のxの点を作成 x_1_grid, x_2_grid = np.meshgrid(x_1_vals, x_2_vals) # 計算用のxの点を作成 x_points = np.stack([x_1_grid.flatten(), x_2_grid.flatten()], axis=1) x_shape = x_1_grid.shape print(x_points.round(3))
[[ 0. 0.513]
[ 0.12 0.513]
[ 0.24 0.513]
...
[11.76 19.487]
[11.88 19.487]
[12. 19.487]]
等高線の設定用の配列を作成します。
# z軸(確率密度)の最小値・最大値を設定 z_min = 0.0 sigma_12_max = np.max([np.abs(sigma_12_vals.min()), np.abs(sigma_12_vals.max())]) z_max = multivariate_t.pdf(x=x_points, loc=mu_d, shape=np.array([[sigma_11, sigma_12_max], [sigma_12_max, sigma_22]]), df=nu).max() z_max = np.ceil(z_max * 10.0) / 10.0 # 等高線を引く値を設定 z_levels = np.linspace(start=z_min, stop=z_max, num=11) print(z_levels)
[0. 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 ]
等高線図のアニメーションを作成します。
# 図を初期化 fig = plt.figure(figsize=(12, 9), facecolor='white') # 図の設定 fig.suptitle("Bivariate Student's t-Distribution", fontsize=20) # 全体のタイトル # カラーバーを表示 tmp = plt.contourf(x_1_grid, x_2_grid, np.zeros(x_shape), vmin=z_min, vmax=z_max, levels=z_levels) # カラーバー用のダミー fig.colorbar(tmp, label='density') # 作図処理を関数として定義 def update(i): # 前フレームのグラフを初期化 plt.cla() # i番目のスケール行列を取得 sigma_12 = sigma_12_vals[i] sigma_dd = np.array([[sigma_11, sigma_12], [sigma_12, sigma_22]]) # 多次元t分布の確率密度を計算 dens = multivariate_t.pdf(x=x_points, loc=mu_d, shape=sigma_dd, df=nu) # パラメータラベル用の文字列を作成 param_text = '$\\nu=' + str(nu) param_text += ', \mu=(' + ', '.join([str(mu) for mu in mu_d]) + ')' param_text += ', \Sigma=' + str([list(sigma_d) for sigma_d in sigma_dd]) + '$' # 2次元t分布を作図:(等高線図) #plt.contour(x_1_grid, x_2_grid, dens.reshape(x_dims), # vmin=z_min, vmax=z_max, levels=z_levels) # 等高線 plt.contourf(x_1_grid, x_2_grid, dens.reshape(x_dims), vmin=z_min, vmax=z_max, levels=z_levels) # 塗りつぶし等高線 plt.xlabel('$x_1$') # x軸ラベル plt.ylabel('$x_2$') # y軸ラベル plt.title(param_text, loc='left') # タイトル plt.grid() # グリッド線 # gif画像を作成 anime_dens = FuncAnimation(fig, update, frames=frame_num, interval=100) # gif画像を保存 anime_dens.save('Bivariate_t_dens_sigma12_cnf.gif')
ヒートマップのアニメーションを作成します。
# 図を初期化 fig = plt.figure(figsize=(12, 9), facecolor='white') # 図の設定 fig.suptitle("Bivariate Student's t-Distribution", fontsize=20) # 全体のタイトル # カラーバーを表示 tmp = plt.pcolor(x_1_grid, x_2_grid, np.zeros(x_shape), vmin=z_min, vmax=z_max) # カラーバー用のダミー fig.colorbar(tmp, label='density') # 作図処理を関数として定義 def update(i): # 前フレームのグラフを初期化 plt.cla() # i番目のスケール行列を取得 sigma_12 = sigma_12_vals[i] sigma_dd = np.array([[sigma_11, sigma_12], [sigma_12, sigma_22]]) # 多次元t分布の確率密度を計算 dens = multivariate_t.pdf(x=x_points, loc=mu_d, shape=sigma_dd, df=nu) # パラメータラベル用の文字列を作成 param_text = '$\\nu=' + str(nu) param_text += ', \mu=(' + ', '.join([str(mu) for mu in mu_d]) + ')' param_text += ', \Sigma=' + str([list(sigma_d) for sigma_d in sigma_dd]) + '$' # 2次元t分布を作図:(ヒートマップ) plt.pcolor(x_1_grid, x_2_grid, dens.reshape(x_dims), vmin=z_min, vmax=z_max) # 塗りつぶし等高線 plt.xlabel('$x_1$') # x軸ラベル plt.ylabel('$x_2$') # y軸ラベル plt.title(param_text, loc='left') # タイトル plt.grid() # グリッド線 # gif画像を作成 anime_dens = FuncAnimation(fig, update, frames=frame_num, interval=100) # gif画像を保存 anime_dens.save('Bivariate_t_dens_sigma12_pcl.gif')
曲面図のアニメーションを作成します。
# 図を初期化 fig = plt.figure(figsize=(9, 9), facecolor='white') # 図の設定 ax = fig.add_subplot(projection='3d') # 3D用の設定 fig.suptitle("Bivariate Student's t-Distribution", fontsize=20) # 全体のタイトル # 作図処理を関数として定義 def update(i): # 前フレームのグラフを初期化 plt.cla() # i番目のスケール行列を取得 sigma_12 = sigma_12_vals[i] sigma_dd = np.array([[sigma_11, sigma_12], [sigma_12, sigma_22]]) # 多次元t分布の確率密度を計算 dens = multivariate_t.pdf(x=x_points, loc=mu_d, shape=sigma_dd, df=nu) # パラメータラベル用の文字列を作成 param_text = '$\\nu=' + str(nu) param_text += ', \mu=(' + ', '.join([str(mu) for mu in mu_d]) + ')' param_text += ', \Sigma=' + str([list(sigma_d) for sigma_d in sigma_dd]) + '$' # 2次元t分布を作図:(曲面図) ax.plot_surface(x_1_grid, x_2_grid, dens.reshape(x_dims), cmap='jet', vmin=z_min, vmax=z_max, alpha=0.8) # 曲面図 ax.contour(x_1_grid, x_2_grid, dens.reshape(x_dims), cmap='jet', vmin=z_min, vmax=z_max, levels=z_levels, offset=0.0) # 等高線図 ax.set_xlabel('$x_1$') # x軸ラベル ax.set_ylabel('$x_2$') # y軸ラベル ax.set_zlabel('density') # z軸ラベル ax.set_title(param_text, loc='left') # タイトル ax.set_zlim(zmin=z_min, zmax=z_max) # z軸の表示範囲 # gif画像を作成 anime_dens = FuncAnimation(fig, update, frames=frame_num, interval=100) # gif画像を保存 anime_dens.save('Bivariate_t_dens_sigma12_srf.gif')
$\sigma_{1,2}$の値に応じて、x軸とy軸の広がる方向が変化するのが分かります。
この記事では、多次元スチューデントのt分布の作図を確認しました。
参考文献
- C.M.ビショップ著,元田 浩・他訳『パターン認識と機械学習 上』,丸善出版,2012年.
おわりに
ほぼ同じコードが何度も出てくることになってしまいました。目が行ったり来たりしないでいいように関数化とかはしないコンセプトなんですが、ここまでくると意思が揺らぎます。とはいうものの、どう関数化すればいいのか思い付いていません。
さて、2022年12月1日は、元Juice=Juiceの宮本佳林さんの24歳のお誕生日です。
今年はソロ活動も本格化して、来年も佳き日々でありますように。そして、私がライブを観に行けますように。
【次の内容】
つづかない
