からっぽのしょこ

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

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

はじめに

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

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

【前の内容】

www.anarchive-beta.com

【他の記事一覧】

www.anarchive-beta.com

【この記事の内容】

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

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

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

# 利用ライブラリ
import numpy as np
from scipy.stats import t
from scipy.special import loggamma, gamma


確率密度の計算

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

標準化t分布

 まずは、標準化されたt分布の計算を行います。

パラメータの設定

 t分布の自由度$\nu$と確率変数の実現値$x$を設定します。

# 自由度を指定
nu = 5

# 確率変数を指定
x = 1.5

 正の整数$\nu$、実数$x$を指定します。設定した値に従う確率密度を計算します。

スクラッチで計算

 定義式から計算します。

# 定義式により確率密度を計算
C = gamma(0.5 * (nu + 1)) / gamma(0.5 * nu) / np.sqrt(np.pi * nu)
term = np.sqrt(1.0 + x**2 / nu)**(nu + 1)
dens = C / term
print(dens)
0.1245173446463551

 t分布は、次の式で定義されます。

$$ \begin{aligned} C_{\mathrm{St}} &= \frac{ \Gamma(\frac{\nu + 1}{2}) }{ \sqrt{\pi \nu} \Gamma(\frac{\nu}{2}) } \\ \mathrm{St}(x | \nu) &= C_{\mathrm{St}} \Bigl( 1 + \frac{x^2}{\nu} \Bigr)^{-\frac{(\nu+1)}{2}} \end{aligned} $$

 ここで、$C_{\mathrm{St}}$はt分布の正規化係数、$\pi$は円周率、$\Gamma(x)$はガンマ関数です。
 円周率はNumPyライブラリのpi、ガンマ関数はSciPyライブラリのgamma()で計算できます。

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

# 対数をとった定義式により確率密度を計算
log_C = loggamma(0.5 * (nu + 1)) - loggamma(0.5 * nu) - 0.5 * np.log(np.pi * nu)
log_term = 0.5 * (nu + 1) * np.log(1.0 + x**2 / nu)
dens = np.exp(log_C - log_term)
print(dens)
0.12451734464635517

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

$$ \begin{aligned} \log C_{\mathrm{St}} &= \log \Gamma \Bigl(\frac{\nu + 1}{2}\Bigr) - \log \Gamma \Bigl(\frac{\nu}{2}\Bigr) - \frac{\log (\pi \nu)}{2} \\ \log \mathrm{St}(x | \nu) &= \log C_{\mathrm{St}} - \frac{(\nu + 1)}{2} \log \left( 1 + \frac{x^2}{\nu} \right) \end{aligned} $$

 対数をとったガンマ関数は、loggamma()で計算できます。引数の値が大きいとgamma()の計算結果が発散してしまいます。その場合でも、loggamma()で計算できます。
 計算結果の指数をとると確率密度が得られます。

$$ \mathrm{St}(x | \nu) = \exp \Bigl( \log \mathrm{St}(x | \nu) \Bigr) $$

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

 次は、SciPyライブラリのクラスを使って確率密度を計算します。

t分布のクラスによる計算

 t分布のクラスtの確率密度メソッドpdf()で計算します。

# t分布の関数により確率密度を計算
dens = t.pdf(x=x, df=nu)
print(dens)
0.12451734464635512

 確率変数の引数xx、自由度の引数dfnuを指定します。

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

# t分布の対数をとった関数により確率密度を計算
log_dens = t.logpdf(x=x, df=nu)
dens = np.exp(log_dens)
print(dens)
0.12451734464635517

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

 他のパラメータのデフォルト値は$\mu = 0, \sigma = 1$です。

# t分布の関数により確率密度を計算
dens = t.pdf(x=x, df=nu, loc=0.0, scale=1.0)
print(dens)
0.12451734464635512

 詳しくは次で確認します。

一般化t分布:スケールパラメータを使用

 次は、スケールパラメータを用いた一般化されたt分布の計算を行います。

パラメータの設定

 t分布のパラメータ$\nu, \mu, \sigma$と確率変数の実現値$x$を設定します。

# 形状パラメータ(自由度)を指定
nu = 5

# 位置パラメータを指定
mu = 2.0

# スケールパラメータを指定
sigma = 0.5

# 確率変数を指定
x = 1.5

 形状パラメータ(自由度)(正の整数)$\nu$、位置パラメータ(実数)$\mu$、スケールパラメータ(正の実数)$\sigma$、実数$x$を指定します。

スクラッチで計算

 定義式から計算します。

# 定義式により確率密度を計算
C = gamma(0.5 * (nu + 1)) / gamma(0.5 * nu)
C /= np.sqrt(np.pi * nu) * sigma
term = np.sqrt(1.0 + ((x - mu) / sigma)**2 / nu)**(nu + 1)
dens = C / term
print(dens)
0.43935959470196134

 $\sigma$を用いたt分布は、次の式で定義されます。

$$ \begin{aligned} C_{\mathrm{St}} &= \frac{ \Gamma(\frac{\nu + 1}{2}) }{ \Gamma(\frac{\nu}{2}) } \frac{1}{\sqrt{\pi \nu} \sigma} \\ \mathrm{St}(x | \nu, \mu, \sigma) &= C_{\mathrm{St}} \left\{ 1 + \frac{1}{\nu} \bigl( \frac{x - \mu}{\sigma} \Bigr)^2 \right\}^{-\frac{(\nu+1)}{2}} \end{aligned} $$

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

# 対数をとった定義式により確率密度を計算
log_C = loggamma(0.5 * (nu + 1)) - loggamma(0.5 * nu)
log_C -= 0.5 * np.log(np.pi * nu) + np.log(sigma)
log_term = 0.5 * (nu + 1) * np.log(1.0 + ((x - mu) / sigma)**2 / nu)
dens = np.exp(log_C - log_term)
print(dens)
0.4393595947019612

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

$$ \begin{aligned} \log C_{\mathrm{St}} &= \log \Gamma \Bigl(\frac{\nu + 1}{2}\Bigr) - \log \Gamma \Bigl(\frac{\nu}{2}\Bigr) - \frac{\log (\pi \nu)}{2} - \log \sigma \\ \log \mathrm{St}(x | \nu, \mu, \sigma) &= \log C_{\mathrm{St}} - \frac{(\nu + 1)}{2} \log \left\{ 1 + \frac{1}{\nu} \Bigl( \frac{x - \mu}{\sigma} \Bigr)^2 \right\} \end{aligned} $$

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

t分布のクラスによる計算

 t分布のクラスtの確率密度メソッドpdf()で計算します。

# t分布の関数により確率密度を計算
dens = t.pdf(x=x, df=nu, loc=mu, scale=sigma)
print(dens)
0.4393595947019612

 確率変数の引数xx、自由度の引数dfnu、位置パラメータの引数locmu、スケールパラメータの引数scalesigmaを指定します。

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

# t分布の対数をとった関数により確率密度を計算
log_dens = t.logpdf(x=x, df=nu, loc=mu, scale=sigma)
dens = np.exp(log_dens)
print(dens)
0.43935959470196123

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

 $\mu, \sigma$の引数を使わずに(標準化t分布により)計算することもできます。

# 標準化t分布の関数により確率密度を計算
y = (x - mu) / sigma # 標準化
dens = t.pdf(x=y, df=nu) / sigma
print(dens)
0.4393595947019612

 次の式で計算できます。

$$ \begin{aligned} y &= \frac{x - \mu}{\sigma} \\ \mathrm{St}(x | \nu, \mu, \sigma) &= \frac{1}{\sigma} \mathrm{St}(y | \nu) \end{aligned} $$

 $x$を$\mu, \sigma$で標準化して$y$とします。自由度$\nu$のときの$y$の確率密度(分布)$\mathrm{St}(y | \nu)$を$\sigma$で割ると、$\nu, \mu, \sigma$のときの$x$の確率密度(分布)$\mathrm{St}(x | \nu, \mu, \sigma)$が得られます。

一般化t分布:逆スケールパラメータを使用

 続いて、逆スケールパラメータを用いた一般化されたt分布の計算を行います。

パラメータの設定

 t分布のパラメータ$\nu, \mu, \lambda$と確率変数の実現値$x$を設定します。

# 形状パラメータ(自由度)を指定
nu = 5

# 位置パラメータを指定
mu = 2.0

# 逆スケールパラメータを指定
lmd = 4.0

# 確率変数を指定
x = 1.5

 スケールパラメータ$\sigma$の代わりに、逆スケールパラメータ(正の実数)$\lambda$を指定します。

 あるいは、$\sigma$を指定して$\lambda$を計算します。

# 逆スケールパラメータを計算
lmd = 1.0 / sigma**2
print(lmd)
4.0

 スケールパラメータ$\sigma$と逆スケールパラメータ$\lambda$は、$\lambda = \frac{1}{\sigma^2}$、$\sigma = \frac{1}{\sqrt{\lambda}}$の関係です。

スクラッチで計算

 定義式から計算します。

# 定義式により確率密度を計算
C = gamma(0.5 * (nu + 1)) / gamma(0.5 * nu)
C *= np.sqrt(lmd / np.pi / nu)
term = np.sqrt(1.0 + lmd / nu * (x - mu)**2)**(nu + 1)
dens = C / term
print(dens)
0.43935959470196134

 $\lambda$を用いたt分布は、次の式で定義されます。

$$ \begin{aligned} C_{\mathrm{St}} &= \frac{ \Gamma(\frac{\nu + 1}{2}) }{ \Gamma(\frac{\nu}{2}) } \Bigl( \frac{\lambda}{\pi \nu} \Bigr)^{\frac{1}{2}} \\ \mathrm{St}(x | \nu, \mu, \lambda) &= C_{\mathrm{St}} \left\{ 1 + \frac{\lambda (x - \mu)^2}{\nu} \right\}^{-\frac{(\nu+1)}{2}} \end{aligned} $$

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

# 対数をとった定義式により確率密度を計算
log_C = loggamma(0.5 * (nu + 1)) - loggamma(0.5 * nu)
log_C -= 0.5 * (np.log(np.pi * nu) - np.log(lmd))
log_term = 0.5 * (nu + 1) * np.log(1.0 + lmd / nu * (x - mu)**2)
dens = np.exp(log_C - log_term)
print(dens)
0.4393595947019612

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

$$ \begin{aligned} \log C_{\mathrm{St}} &= \log \Gamma \Bigl(\frac{\nu + 1}{2}\Bigr) - \log \Gamma \Bigl(\frac{\nu}{2}\Bigr) - \frac{1}{2} \Bigl\{ \log (\pi \nu) - \log \lambda \Bigr\} \\ \log \mathrm{St}(x | \nu, \mu, \lambda) &= \log C_{\mathrm{St}} - \frac{(\nu + 1)}{2} \log \left\{ 1 + \frac{\lambda (x - \mu)^2}{\nu} \right\} \end{aligned} $$

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

t分布のクラスによる計算

 pdf()またはlogpdf()で計算します。

# t分布の関数により確率密度を計算
dens = t.pdf(x=x, df=nu, loc=mu, scale=1.0/np.sqrt(lmd))
print(dens)

# t分布の対数をとった関数により確率密度を計算
log_dens = t.logpdf(x=x, df=nu, loc=mu, scale=1.0/np.sqrt(lmd))
dens = np.exp(log_dens)
print(dens)
0.4393595947019612
0.43935959470196123

 スケールパラメータの引数scalelmdの平方根の逆数を指定します。

 $\mu, \sigma$の引数を使わずに計算することもできます。

# 標準化t分布の関数により確率密度を計算
y = (x - mu) * np.sqrt(lmd) # 標準化
dens = t.pdf(x=y, df=nu) * np.sqrt(lmd)
print(dens)
0.4393595947019612

 次の式で計算できます。

$$ \begin{aligned} y &= (x - \mu) \sqrt{\lambda} \\ \mathrm{St}(x | \nu, \mu, \lambda) &= \sqrt{\lambda} \mathrm{St}(y | \nu) \end{aligned} $$

 $x$を$\mu, \lambda$で標準化して$y$とします。自由度$\nu$のときの$y$の確率密度(分布)$\mathrm{St}(y | \nu)$に$\lambda$の平方根を掛けると、$\nu, \mu, \lambda$のときの$x$の確率密度(分布)$\mathrm{St}(x | \nu, \mu, \lambda)$が得られます。
 「スケールパラメータを使用」のときの計算式について、$\frac{1}{\sigma} = \sqrt{\lambda}$で置き換えた式です。

統計量の計算

 次は、スチューデントのt分布の期待値・分散・最頻値を計算します。詳しくはいつか書きます。

スクラッチで計算

 期待値を計算します。

# 計算式により期待値を計算:(nu > 1)
E_x = mu
print(E_x)
2.0

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

 分散を計算します。

# 計算式により分散を計算:(nu > 2)
V_x = sigma**2 * nu / (nu - 2.0)
print(V_x)

# 計算式により分散を計算:(nu > 2)
V_x = nu / (nu - 2.0) / lmd
print(V_x)
0.4166666666666667
0.4166666666666667

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

$$ \mathbb{V}[x] = \sigma^2 \frac{\nu}{\nu - 2} = \frac{1}{\lambda} \frac{\nu}{\nu - 2} $$

 $\sigma^2 = \frac{1}{\lambda}$が成り立つのが分かります。

 最頻値を計算します。

# 計算式により最頻値を計算
mode_x = mu
print(mode_x)
2.0

 t分布の最頻値は$\mu$です。こちらは、$\nu$に関わらず計算できます。

 次は、SciPyライブラリのクラスを使って統計量を計算します。

t分布のクラスによる計算

 期待値メソッドmean()で計算します。

# t分布の関数により期待値を計算:(nu > 1)
E_x = t.mean(df=nu, loc=mu, scale=sigma)
print(E_x)
2.0

 確率密度メソッドと同様に引数を指定します。

 分散メソッドvar()で計算します。

# t分布の関数により分散を計算:(nu > 2)
V_x = t.var(df=nu, loc=mu, scale=sigma)
print(V_x)
0.4166666666666667

 こちらも同様です。

 この記事では、1次元スチューデントのt分布の計算を確認しました。次は、グラフを作成します。

参考文献

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

おわりに

 見て見ぬふりしてたスケールパラメータについて解決したので満足しました。

【次の内容】

www.anarchive-beta.com