はじめに
機械学習で登場する確率分布について色々な角度から理解したいシリーズです。
この記事では、R言語で1次元(一変量)スチューデントのt分布に関する計算をします。
【前の内容】
【他の記事一覧】
【この記事の内容】
1次元スチューデントのt分布の計算
1次元スチューデントのt分布(Student's t-Distribution)の確率密度と統計量を計算します。t分布については「1次元スチューデントのt分布の定義式の確認 - からっぽのしょこ」を参照してください。
利用するパッケージを読み込みます。
# 利用パッケージ library(LaplacesDemon)
この記事では、パッケージ名::関数名()の記法を使うので、パッケージを読み込む必要はありません。
一般化したt分布の計算にLaplacesDemonパッケージを利用します。
確率密度の計算
スチューデントのt分布に従う確率密度を計算する方法をいくつか確認します。
標準化t分布
まずは、標準化されたt分布の計算を行います。
パラメータの設定
t分布の自由度$\nu$と確率変数の実現値$x$を設定します。
# 自由度を指定 nu <- 5 # 確率変数を指定 x <- 1.5
正の整数$\nu$、実数$x$を指定します。設定した値に従う確率密度を計算します。
スクラッチで計算
定義式から計算します。
# 定義式により確率密度を計算 C <- gamma(0.5 * (nu + 1)) / gamma(0.5 * nu) / sqrt(pi * nu) term <- sqrt(1 + x^2 / nu)^(nu + 1) dens <- C / term dens
## [1] 0.1245173
t分布は、次の式で定義されます。
ここで、$C_{\mathrm{St}}$はt分布の正規化係数、$\pi$は円周率、$\Gamma(x)$はガンマ関数です。
円周率はpi、ガンマ関数はgamma()で計算できます。
対数をとった定義式から計算します。
# 対数をとった定義式により確率密度を計算 log_C <- lgamma(0.5 * (nu + 1)) - lgamma(0.5 * nu) - 0.5 * log(pi * nu) log_term <- 0.5 * (nu + 1) * log(1 + x^2 / nu) log_dens <- log_C - log_term dens <- exp(log_dens) dens; log_dens
## [1] 0.1245173 ## [1] -2.08331
対数をとった定義式を計算します。
対数をとったガンマ関数は、lgamma()で計算できます。引数の値が大きいとgamma()の計算結果が発散してしまいます。その場合でも、lgamma()で計算できます。
計算結果の指数をとると確率密度が得られます。
指数と対数の性質より$\exp(\log x) = x$です。
次は、関数を使って確率密度を計算します。
t分布の関数による計算
t分布の確率密度関数dt()で計算します。
# 標準化t分布の関数により確率密度を計算 dens <- dt(x = x, df = nu) dens
## [1] 0.1245173
確率変数の引数xにx、自由度の引数dfにnuを指定します。
log = TRUEを指定すると対数をとった確率密度を返します。
# 標準化t分布の対数をとった関数により確率密度を計算 log_dens <- dt(x = x, df = nu, log = TRUE) dens <- exp(log_dens) dens; log_dens
## [1] 0.1245173 ## [1] -2.08331
計算結果の指数をとると確率密度が得られます。
一般化t分布：スケールパラメータを使用
次は、スケールパラメータを用いた一般化されたt分布の計算を行います。
パラメータの設定
t分布のパラメータ$\nu, \mu, \sigma$と確率変数の実現値$x$を設定します。
# 形状パラメータ(自由度)を指定 nu <- 5 # 位置パラメータを指定 mu <- 2 # スケールパラメータを指定 sigma <- 0.5 # 確率変数を指定 x <- 1.5
形状パラメータ(自由度)(正の整数)$\nu$、位置パラメータ(実数)$\mu$、スケールパラメータ(正の実数)$\sigma$、実数$x$を指定します。
スクラッチで計算
定義式から計算します。
# 対数をとった定義式により確率密度を計算 C <- gamma(0.5 * (nu + 1)) / gamma(0.5 * nu) / sqrt(pi * nu) / sigma term <- sqrt(1 + ((x - mu) / sigma)^2 / nu)^(nu + 1) dens <- C / term dens
## [1] 0.4393596
$\sigma$を用いたt分布は、次の式で定義されます。
対数をとった定義式から計算します。
# 対数をとった定義式により確率密度を計算 log_C <- lgamma(0.5 * (nu + 1)) - lgamma(0.5 * nu) - 0.5 * log(pi * nu) - log(sigma) log_term <- 0.5 * (nu + 1) * log(1 + ((x - mu) / sigma)^2 / nu) log_dens <- log_C - log_term dens <- exp(log_dens) dens; log_dens
## [1] 0.4393596 ## [1] -0.8224371
対数をとった定義式を計算します。
計算結果の指数をとると確率密度が得られます。
t分布のクラスによる計算
LaplacesDemonパッケージの一般化t分布の確率密度関数dst()で計算します。
# 一般化t分布の関数により確率密度を計算 dens <- LaplacesDemon::dst(x = x, mu = mu, sigma = sigma, nu = nu) dens
## [1] 0.4393596
確率変数の引数xにx、位置パラメータの引数muにmu、スケールパラメータの引数sigmaにsigma、形状パラメータ(自由度)の引数nuにnuを指定します。
log = TRUEを指定すると対数をとった確率密度を計算します。
# 一般化t分布の対数をとった関数により確率密度を計算 log_dens <- LaplacesDemon::dst(x = x, mu = mu, sigma = sigma, nu = nu, log = TRUE) dens <- exp(log_dens) dens; log_dens
## [1] 0.4393596 ## [1] -0.8224371
計算結果の指数をとると確率密度が得られます。
標準化t分布の関数dt()でも計算できます。
# 標準化t分布の関数により確率密度を計算 y <- (x - mu) / sigma # 標準化 dens <- dt(x = y, df = nu) / sigma dens
## [1] 0.4393596
次の式で計算できます。
$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 # 逆スケールパラメータを指定 lambda <- 4 # 確率変数を指定 x <- 1.5
スケールパラメータ$\sigma$の代わりに、逆スケールパラメータ(正の実数)$\lambda$を指定します。
あるいは、$\sigma$を指定して$\lambda$を計算します。
# 逆スケールパラメータを計算 lambda <- 1 / sigma^2 lambda
## [1] 4
スケールパラメータ$\sigma$と逆スケールパラメータ$\lambda$は、$\lambda = \frac{1}{\sigma^2}$、$\sigma = \frac{1}{\sqrt{\lambda}}$の関係です。
スクラッチで計算
定義式から計算します。
# 定義式により確率密度を計算 C <- gamma(0.5 * (nu + 1)) / gamma(0.5 * nu) * sqrt(lambda / pi / nu) term <- sqrt(1 + lambda * (x - mu)^2 / nu)^(nu + 1) dens <- C / term dens
## [1] 0.4393596
$\lambda$を用いたt分布は、次の式で定義されます。
対数をとった定義式から計算します。
# 対数をとった定義式により確率密度を計算 log_C <- lgamma(0.5 * (nu + 1)) - lgamma(0.5 * nu) - 0.5 * (log(pi * nu) - log(lambda)) log_term <- 0.5 * (nu + 1) * log(1 + lambda / nu * (x - mu)^2) log_dens <- log_C - log_term dens <- exp(log_dens) dens; log_dens
## [1] 0.4393596 ## [1] -0.8224371
対数をとった定義式を計算します。
計算結果の指数をとると確率密度が得られます。
t分布のクラスによる計算
LaplacesDemonパッケージのdst()で計算します。
# 一般化t分布の関数により確率密度を計算 dens <- LaplacesDemon::dst(x = x, mu = mu, sigma = 1/sqrt(lambda), nu = nu) dens # 一般化t分布の対数をとった関数により確率密度を計算 log_dens <- LaplacesDemon::dst(x = x, mu = mu, sigma = 1/sqrt(lambda), nu = nu, log = TRUE) dens <- exp(log_dens) dens; log_dens
## [1] 0.4393596 ## [1] 0.4393596 ## [1] -0.8224371
スケールパラメータの引数sigmaにlambdaの平方根の逆数を指定します。
$\mu, \sigma$の引数を使わずに計算することもできます。
# 標準化t分布の関数により確率密度を計算 y <- (x - mu) * sqrt(lambda) # 標準化 dens <- dt(x = y, df = nu) * sqrt(lambda) dens
## [1] 0.4393596
次の式で計算できます。
$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 E_x
## [1] 2
t分布の期待値は$\mu$です。ただし、$\nu > 1$の場合に定義されます。
分散を計算します。
# sigmaを使って分散を計算:(nu > 2) V_x <- sigma^2 * nu / (nu - 2) V_x # lambdaを使って分散を計算:(nu > 2) V_x <- nu / (nu - 2) / lambda V_x
## [1] 0.4166667 ## [1] 0.4166667
t分布の分散は、$\sigma$または$\lambda$を使って、次の式で計算できます。ただし、$\nu > 2$の場合に定義されます。
$\sigma^2 = \frac{1}{\lambda}$が成り立つのが分かります。
最頻値を計算します。
# 最頻値を計算 mode_x <- mu mode_x
## [1] 2
t分布の最頻値は$\mu$です。こちらは、$\nu$に関わらず計算できます。
この記事では、1次元スチューデントのt分布の計算を確認しました。次は、グラフを作成します。
参考文献
- C.M.ビショップ著,元田 浩・他訳『パターン認識と機械学習 上』,丸善出版,2012年.
おわりに
半分コピペで済むので楽なのですが、楽だと楽しくなくて全然進まないです。
【次の内容】
