はじめに
機械学習で登場する確率分布について色々な角度から理解したいシリーズです。
この記事では、Pythonで1次元(一変量)スチューデントのt分布に関する計算をします。
【前の内容】
【他の記事一覧】
【この記事の内容】
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分布は、次の式で定義されます。
ここで、$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
対数をとった定義式を計算します。
対数をとったガンマ関数は、loggamma()
で計算できます。引数の値が大きいとgamma()
の計算結果が発散してしまいます。その場合でも、loggamma()
で計算できます。
計算結果の指数をとると確率密度が得られます。
指数と対数の性質より$\exp(\log x) = x$です。
次は、SciPy
ライブラリのクラスを使って確率密度を計算します。
t分布のクラスによる計算
t分布のクラスt
の確率密度メソッドpdf()
で計算します。
# t分布の関数により確率密度を計算 dens = t.pdf(x=x, df=nu) print(dens)
0.12451734464635512
確率変数の引数x
にx
、自由度の引数df
にnu
を指定します。
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分布は、次の式で定義されます。
対数をとった定義式から計算します。
# 対数をとった定義式により確率密度を計算 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
対数をとった定義式を計算します。
計算結果の指数をとると確率密度が得られます。
t分布のクラスによる計算
t分布のクラスt
の確率密度メソッドpdf()
で計算します。
# t分布の関数により確率密度を計算 dens = t.pdf(x=x, df=nu, loc=mu, scale=sigma) print(dens)
0.4393595947019612
確率変数の引数x
にx
、自由度の引数df
にnu
、位置パラメータの引数loc
にmu
、スケールパラメータの引数scale
にsigma
を指定します。
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
次の式で計算できます。
$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分布は、次の式で定義されます。
対数をとった定義式から計算します。
# 対数をとった定義式により確率密度を計算 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
対数をとった定義式を計算します。
計算結果の指数をとると確率密度が得られます。
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
スケールパラメータの引数scale
にlmd
の平方根の逆数を指定します。
$\mu, \sigma$の引数を使わずに計算することもできます。
# 標準化t分布の関数により確率密度を計算 y = (x - mu) * np.sqrt(lmd) # 標準化 dens = t.pdf(x=y, df=nu) * np.sqrt(lmd) print(dens)
0.4393595947019612
次の式で計算できます。
$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$の場合に定義されます。
$\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年.
おわりに
見て見ぬふりしてたスケールパラメータについて解決したので満足しました。
【次の内容】