からっぽのしょこ

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

トップ > 確率分布 > 【Python】多次元スチューデントのt分布の計算

【Python】多次元スチューデントのt分布の計算

確率分布 確率分布-スチューデントのt分布 Python

はじめに

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

　この記事では、Pythonで多次元(多変量)スチューデントのt分布に関する計算をします。

【前の内容】

www.anarchive-beta.com

【他の記事一覧】

www.anarchive-beta.com

【この記事の内容】

多次元スチューデントのt分布の計算

　多次元スチューデントのt分布(Multivariate Student's t Distribution)の確率密度を計算します。多次元t分布については「多次元スチューデントのt分布の定義式 - からっぽのしょこ」を参照してください。

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

# ライブラリを読み込み
import numpy as np
from scipy.special import gamma, loggamma
from scipy.stats import multivariate_t


確率密度の計算

　多次元スチューデントのt分布に従う確率密度を計算する方法をいくつか確認します。

スケール行列を使用

　まずは、スケール行列を用いてt分布を計算します。

パラメータの設定

　t分布の次元数$D$とパラメータ$\boldsymbol{\mu}, \boldsymbol{\Sigma}$、確率変数の値$\mathbf{x}$を設定します。

# 次元数を指定
D = 3

# 自由度を指定
nu = 5

# 位置ベクトルを指定
mu_d = np.array([10.0, -6.0, 1.5])

# スケール行列を指定
sigma_dd = np.array(
    [[4.0, 1.8, -0.1], 
     [1.8, 9.0, 2.4], 
     [-0.1, 2.4, 1.0]]
)

# 確率変数の値を指定
x_d = np.array([11.5, -5.0, 0.0])

　自由度(形状パラメータ)$\nu$、位置ベクトル$\boldsymbol{\mu}$、スケール行列$\boldsymbol{\Sigma}$、確率変数の値$\mathbf{x}$を指定します。$\boldsymbol{\mu}, \mathbf{x}$は$D$次元ベクトル、$\boldsymbol{\Sigma}$は$D \times D$の行列です。

$$ \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} ,\ \mathbf{x} = \begin{pmatrix} x_1 \\ x_2 \\ \vdots \\ x_D \end{pmatrix} $$

　$x_d$は実数をとり、$\nu$は正の実数、$\mu_d$は実数、$\sigma_{d,d}$は正の実数、$\sigma_{i,j} = \sigma_{j,i}\ (i \neq j)$は実数、また$\boldsymbol{\Sigma}$は正定値行列を満たす必要があります。設定した値に従う確率密度を計算します。

スクラッチで計算

　定義式から計算します。

# 定義式により確率密度を計算
C_St = gamma(0.5 * (nu + D)) / gamma(0.5 * nu)
C_St /= np.sqrt(np.pi * nu)**D * np.sqrt(np.linalg.det(sigma_dd))
x_tilde_d1 = (x_d - mu_d).reshape((D, 1))
tmp_term = x_tilde_d1.T.dot(np.linalg.inv(sigma_dd)).dot(x_tilde_d1).item()
dens = C_St / np.sqrt(1.0 + tmp_term / nu)**(nu + D)
print(dens)

0.0003309268920926894

　t分布は、$\boldsymbol{\Sigma}$を用いて次の式で定義されます。

$$ \begin{aligned} C_{\mathrm{St}}(\nu, \boldsymbol{\Sigma}) &= \frac{ \Gamma(\frac{\nu + D}{2}) }{ \Gamma(\frac{\nu}{2}) } \frac{ 1 }{ (\pi \nu)^{\frac{D}{2}} |\boldsymbol{\Sigma}|^{\frac{1}{2}} } \\ \mathrm{St}(\mathbf{x} | \nu, \boldsymbol{\mu}, \boldsymbol{\Sigma}) &= C_{\mathrm{St}}(\nu, \boldsymbol{\Sigma}) \left\{ 1 + \frac{1}{\nu} (\mathbf{x} - \boldsymbol{\mu})^{\top} \boldsymbol{\Sigma}^{-1} (\mathbf{x} - \boldsymbol{\mu}) \right\}^{-\frac{\nu+D}{2}} \end{aligned} $$

　ここで、$C_{\mathrm{St}}$はt分布の正規化係数、$\pi$は円周率、$\mathbf{A}^{\top}$は転置行列、$\mathbf{A}^{-1}$は逆行列、$|\mathbf{A}|$は行列式、2分の1乗は平方根$\sqrt{a} = a^{\frac{1}{2}}$です。
　円周率はnp.pi、転置は.Tメソッド、逆行列np.linalg.inv()、行列式はnp.linalg.det()、行列の積はnp.dot()で計算できます。

　対数をとった定義式から計算します。

# 対数をとった定義式により確率密度を計算
log_C_St = loggamma(0.5 * (nu + D)) - loggamma(0.5 * nu)
log_C_St -= D * 0.5 * np.log(np.pi * nu) + 0.5 * np.log(np.linalg.det(sigma_dd))
x_tilde_d1 = (x_d - mu_d).reshape((D, 1))
tmp_term = x_tilde_d1.T.dot(np.linalg.inv(sigma_dd)).dot(x_tilde_d1).item()
log_dens = log_C_St - (nu + D) * 0.5 * np.log(1.0 + tmp_term / nu)
dens = np.exp(log_dens)
print(dens, log_dens)

0.0003309268920926896 -8.013613076793145

　対数をとった定義式を計算します。

$$ \begin{aligned} \log C_{\mathrm{St}}(\nu, \boldsymbol{\Sigma}) &= \log \Gamma \Bigl(\frac{\nu + D}{2}\Bigr) - \log \Gamma \Bigl(\frac{\nu}{2}\Bigr) - \frac{D}{2} \log (\pi \nu) - \frac{1}{2} \log |\boldsymbol{\Sigma}| \\ \log \mathrm{St}(\mathbf{x} | \nu, \boldsymbol{\mu}, \boldsymbol{\Sigma}) &= \log C_{\mathrm{St}}(\nu, \boldsymbol{\Sigma}) - \frac{\nu+D}{2} \log \Bigl\{ 1 + \frac{1}{\nu} (\mathbf{x} - \boldsymbol{\mu})^{\top} \boldsymbol{\Sigma}^{-1} (\mathbf{x} - \boldsymbol{\mu}) \Bigr\} \end{aligned} $$

　計算結果の指数をとると確率密度が得られます。

$$ \mathrm{St}(\mathbf{x} | \nu, \boldsymbol{\mu}, \boldsymbol{\Sigma}) = \exp \Bigr( \log \mathrm{St}(\mathbf{x} | \nu, \boldsymbol{\mu}, \boldsymbol{\Sigma}) \Bigr) $$

　指数と対数の性質より$\exp(\log x) = x$です。

　次は、モジュールを使って確率密度を計算します。

クラスで計算

　SciPyライブラリのstatsモジュールの多次元t分布のクラスmultivariate_tの確率密度メソッドpdf()で計算します。

# 関数により確率密度を計算
dens = multivariate_t.pdf(x=x_d, loc=mu_d, shape=sigma_dd, df=nu)
print(dens)

0.00033092689209273426

　確率変数の引数xにx_d、位置ベクトルの引数locにmu_d、スケール行列の引数shapeにsigma_dd、自由度の引数dfにnuを指定します。

　logpdf()だと対数をとった確率密度を返します。

# 対数をとった関数により確率密度を計算
log_dens = multivariate_t.logpdf(x=x_d, loc=mu_d, shape=sigma_dd, df=nu)
print(dens, log_dens)

0.00033092689209273426 -8.01361307679301

　計算結果の指数をとると確率密度が得られます。

逆スケール行列を使用

　次は、逆スケール行列を用いてt分布を計算します。

パラメータの設定

　スケール行列$\boldsymbol{\Sigma}$の代わりに、逆スケール行列$\boldsymbol{\Lambda}$を指定します。

# 逆スケール行列を指定
lambda_dd = np.array(
    [[4.0, 1.8, -0.1], 
     [1.8, 9.0, 2.4], 
     [-0.1, 2.4, 1.0]]
)


　あるいは、$\boldsymbol{\Sigma}$を指定して$\boldsymbol{\Lambda}$を計算します。

# 逆スケール行列を計算
lambda_dd = np.linalg.inv(sigma_dd)
print(lambda_dd)

[[ 0.36960986 -0.23271732  0.59548255]
 [-0.23271732  0.45516769 -1.1156742 ]
 [ 0.59548255 -1.1156742   3.73716632]]

　逆スケール行列は、スケール行列の逆行列$\boldsymbol{\Lambda} = \boldsymbol{\Sigma}^{-1}$です。

スクラッチで計算

　定義式から計算します。

# 定義式により確率密度を計算
C_St = gamma(0.5 * (nu + D)) / gamma(0.5 * nu)
C_St *= np.sqrt(np.linalg.det(lambda_dd)) / np.sqrt(np.pi * nu)**D
x_tilde_d1 = (x_d - mu_d).reshape((D, 1))
tmp_term = x_tilde_d1.T.dot(lambda_dd).dot(x_tilde_d1).item()
dens = C_St / np.sqrt(1.0 + tmp_term / nu)**(nu + D)
print(dens)

0.0003309268920926894

　t分布は、$\boldsymbol{\Lambda}$を用いて次の式で定義されます。

$$ \begin{aligned} C_{\mathrm{St}}(\nu, \boldsymbol{\Lambda}) &= \frac{ \Gamma(\frac{\nu + D}{2}) }{ \Gamma(\frac{\nu}{2}) } \frac{ |\boldsymbol{\Lambda}|^{\frac{1}{2}} }{ (\pi \nu)^{\frac{D}{2}} } \\ \mathrm{St}(\mathbf{x} | \nu, \boldsymbol{\mu}, \boldsymbol{\Lambda}) &= C_{\mathrm{St}}(\nu, \boldsymbol{\Lambda}) \left\{ 1 + \frac{1}{\nu} (\mathbf{x} - \boldsymbol{\mu})^{\top} \boldsymbol{\Lambda} (\mathbf{x} - \boldsymbol{\mu}) \right\}^{-\frac{\nu+D}{2}} \end{aligned} $$

　対数をとった定義式から計算します。

# 対数をとった定義式により確率密度を計算
log_C_St = loggamma(0.5 * (nu + D)) - loggamma(0.5 * nu)
log_C_St += 0.5 * np.log(np.linalg.det(lambda_dd)) - D * 0.5 * np.log(np.pi * nu)
x_tilde_d1 = (x_d - mu_d).reshape((D, 1))
tmp_term = x_tilde_d1.T.dot(lambda_dd).dot(x_tilde_d1).item()
log_dens = log_C_St - (nu + D) * 0.5 * np.log(1.0 + tmp_term / nu)
dens = np.exp(log_dens)
print(dens, log_dens)

0.0003309268920926896 -8.013613076793145

　対数をとった定義式を計算します。

$$ \begin{aligned} \log C_{\mathrm{St}}(\nu, \boldsymbol{\Lambda}) &= \log \Gamma \Bigl(\frac{\nu + D}{2}\Bigr) - \log \Gamma \Bigl(\frac{\nu}{2}\Bigr) - \frac{D}{2} \log (\pi \nu) + \frac{1}{2} \log |\boldsymbol{\Lambda}| \\ \log \mathrm{St}(\mathbf{x} | \nu, \boldsymbol{\mu}, \boldsymbol{\Lambda}) &= \log C_{\mathrm{St}}(\nu, \boldsymbol{\Lambda}) - \frac{\nu+D}{2} \log \Bigl\{ 1 + \frac{1}{\nu} (\mathbf{x} - \boldsymbol{\mu})^{\top} \boldsymbol{\Lambda} (\mathbf{x} - \boldsymbol{\mu}) \Bigr\} \end{aligned} $$

　計算結果の指数をとると確率密度が得られます。

クラスで計算

　multivariate_tクラスのpdf()またはlogpdf()で計算します。

# 関数により確率密度を計算
dens = multivariate_t.pdf(x=x_d, loc=mu_d, shape=np.linalg.inv(lambda_dd), df=nu)
print(dens)

# 対数をとった関数により確率密度を計算
log_dens = multivariate_t.logpdf(x=x_d, loc=mu_d, shape=np.linalg.inv(lambda_dd), df=nu)
print(dens, log_dens)

0.00033092689209272954
0.00033092689209272954 -8.013613076793025

　スケール行列の引数shapeにlambda_ddの逆行列を指定します。

標準化t分布による計算

　続いて、確率変数の値を変換しておき、標準化t分布により確率密度を計算します。スケール行列と固有値・固有ベクトルの関係については「2.3.0：分散共分散行列と固有値・固有ベクトルの関係の導出【PRMLのノート】 - からっぽのしょこ」を参照してください。

　スケール行列の固有値と固有ベクトルを計算します。

# 固有値・固有ベクトルを計算
res_eigen = np.linalg.eig(sigma_dd)
lambda_d, u_dd = res_eigen[0], res_eigen[1].T
print(lambda_d)
print(u_dd)

[10.14233952  3.61882683  0.23883365]
[[ 0.26925095  0.93222699  0.24177834]
 [ 0.94887152 -0.21383579 -0.23220055]
 [ 0.16476276 -0.2919368   0.94213913]]

　固有値$\lambda_d$と固有ベクトル$\mathbf{u}_d = (u_{d,1}, \cdots, u_{d, D})$をnp.linalg.eig()で計算します。タプルが出力されるので、0番目の要素の固有値(をまとめた1次元配列)、1番目の要素の固有ベクトル(をまとめた2次元配列)を取り出します。数式での成分と合わせるために転置しておきます。

　確率変数の値を変換します。

# 変数を変換
y_d = np.dot(u_dd, (x_d - mu_d).reshape((D, 1))).flatten() / np.sqrt(lambda_d)
print(y_d)

[ 0.30565977  0.81887966 -2.98339094]

　$D$個の固有値$\lambda_d$を対角成分とする対角行列を$\mathbf{V}$(固有値の逆数$\frac{1}{\lambda_d}$を対角成分とする対角行列を$\mathbf{V}^{-1}$)、固有ベクトル$\mathbf{u}_d^{\top}$を行方向に並べた行列を$\mathbf{U}$として、次の式を計算します。

$$ \begin{aligned} \mathbf{y} &= \mathbf{V}^{-\frac{1}{2}} \mathbf{U} (\mathbf{x} - \boldsymbol{\mu}) \\ y_d &= \frac{1}{\sqrt{\lambda_d}} \mathbf{u}_d^{\top} (\mathbf{x} - \boldsymbol{\mu}) \end{aligned} $$

　$\mathbf{x}$の確率密度を計算します。

# Σを用いて標準化t分布により確率密度を計算
dens = multivariate_t.pdf(x=y_d, loc=np.zeros(D), shape=np.identity(D), df=nu)
dens /= np.sqrt(np.linalg.det(sigma_dd))
print(dens)

0.00033092689209268997

　標準化t分布に従う$\mathbf{y}$の確率密度を計算して、$\boldsymbol{\Sigma}$の行列式の平方根の逆数を掛けます。

$$ \mathrm{St}(\mathbf{x} | \nu, \boldsymbol{\mu}, \boldsymbol{\Sigma}) = \frac{ 1 }{ |\boldsymbol{\Sigma}|^{\frac{1}{2}} } \mathrm{St}(\mathbf{y} | \nu) $$

　標準化t分布に従う$\mathbf{y}$の確率密度を計算して、$\boldsymbol{\Sigma}$の行列式の平方根の逆数を掛けます。
　標準化t分布の計算は、multivariate_t.pdf()の位置ベクトルの引数locに0ベクトル、スケール行列の引数shapeに単位行列を指定して計算できます。

　逆スケール行列を用いる場合は、次のように計算します。

# Λを用いて標準化t分布により確率密度を計算
dens = np.sqrt(np.linalg.det(lambda_dd))
dens *= multivariate_t.pdf(x=y_d, loc=np.zeros(D), shape=np.identity(D), df=nu)
print(dens)

0.0003309268920926901

　標準化t分布に従う$\mathbf{y}$の確率密度を計算して、$\boldsymbol{\Lambda}$の行列式の平方根を掛けます。

$$ \mathrm{St}(\mathbf{x} | \nu, \boldsymbol{\mu}, \boldsymbol{\Lambda}) = |\boldsymbol{\Lambda}|^{\frac{1}{2}} \mathrm{St}(\mathbf{y} | \nu) $$


統計量の計算

　最後に、多次元スチューデントのt分布の期待値・共分散・最頻値を計算します。

　期待値を計算します。

# 期待値を計算:(ν > 1)
E_x_d = mu_d
print(E_x_d)

[10.  -6.   1.5]

　t分布の期待値は$\boldsymbol{\mu}$です。ただし、$\nu > 1$の場合に定義されます。

$$ \mathbb{E}[\mathbf{x}] = \boldsymbol{\mu} \quad (\nu > 1) $$

　共分散を計算します。

# Σを使って共分散を計算:(ν > 2)
cov_x_dd = nu / (nu - 2) * sigma_dd
print(cov_x_dd)

# Λを使って共分散を計算:(ν > 2)
cov_x_dd = nu / (nu - 2) * np.linalg.inv(lambda_dd)
print(cov_x_dd)

[[ 6.66666667  3.         -0.16666667]
 [ 3.         15.          4.        ]
 [-0.16666667  4.          1.66666667]]
[[ 6.66666667  3.         -0.16666667]
 [ 3.         15.          4.        ]
 [-0.16666667  4.          1.66666667]]

　t分布の分散は、$\boldsymbol{\Sigma}$または$\boldsymbol{\Lambda}$を使って、次の式で計算できます。ただし、$\nu > 2$の場合に定義されます。

$$ \mathrm{Cov}[\mathbf{x}] = \frac{\nu}{\nu - 2} \boldsymbol{\Sigma} = \frac{\nu}{\nu - 2} \boldsymbol{\Lambda}^{-1} \quad (\nu > 2) $$

　$\boldsymbol{\Sigma} = \boldsymbol{\Lambda}^{-1}$が成り立つのが分かります。

　最頻値を計算します。

# 最頻値を計算
mode_x_d = mu_d
print(mode_x_d)

[10.  -6.   1.5]

　t分布の最頻値は$\boldsymbol{\mu}$です。

$$ \mathrm{mode}[\mathbf{x}] = \boldsymbol{\mu} $$


　この記事では、多次元t分布の計算を確認しました。次は、グラフを作成します。

参考文献

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

おわりに

　今日からこのブログの開設5年目で、この記事がその1つ目です。まだ続けられる。

【次の内容】

www.anarchive-beta.com