からっぽのしょこ

読んだら書く!書いたら読む!同じ事は二度調べ(たく)ない

【Python】2次元スチューデントのt分布の作図

はじめに

 機械学習で登場する確率分布について色々な角度から理解したいシリーズです。

 この記事では、Pythonで多次元(多変量)スチューデントのt分布のグラフを作成します。

【前の内容】

www.anarchive-beta.com

【他の記事一覧】

www.anarchive-beta.com

【この記事の内容】

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_valsx_2_valsの要素の全ての組み合わせ(格子状の点)をnp.meshgrid()で作成してx_1_grid, x_2_gridとします。
 x_1_grid, x_2_gridnp.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()で計算できます。確率変数の引数xx_points、位置ベクトルの引数locmu_d、スケール行列の引数shapesigma_dd、自由度の引数dfnuを指定します。

 パラメータの値を数式で表示するための文字列を作成します。

# パラメータラベル用の文字列を作成
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() # 描画

2次元スチューデントのt分布の等高線図

 等高線は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() # 描画

2次元スチューデントのt分布のヒートマップ

 ヒートマップは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() # 描画

2次元スチューデントのt分布の3Dグラフ

 曲面図は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')


2次元スチューデントのt分布の自由度と形状の関係

 $\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')


2次元スチューデントのt分布の位置ベクトルと形状の関係

 $\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')


2次元スチューデントのt分布の位置ベクトルと形状の関係

 $\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')


2次元スチューデントのt分布のスケール行列と形状の関係

 $\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')


2次元スチューデントのt分布のスケール行列と形状の関係

 $\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')


2次元スチューデントのt分布のスケール行列と形状の関係

 $\sigma_{1,2}$の値に応じて、x軸とy軸の広がる方向が変化するのが分かります。

 この記事では、多次元スチューデントのt分布の作図を確認しました。

参考文献

  • C.M.ビショップ著,元田 浩・他訳『パターン認識と機械学習 上』,丸善出版,2012年.

おわりに

 ほぼ同じコードが何度も出てくることになってしまいました。目が行ったり来たりしないでいいように関数化とかはしないコンセプトなんですが、ここまでくると意思が揺らぎます。とはいうものの、どう関数化すればいいのか思い付いていません。

 さて、2022年12月1日は、元Juice=Juiceの宮本佳林さんの24歳のお誕生日です。

 今年はソロ活動も本格化して、来年も佳き日々でありますように。そして、私がライブを観に行けますように。

【次の内容】

つづかない