はじめに
機械学習で登場する確率分布について色々な角度から理解したいシリーズです。
この記事では、R言語で1次元(一変量)スチューデントのt分布のグラフを作成します。
【前の内容】
【他の記事一覧】
【この記事の内容】
1次元スチューデントのt分布の作図
1次元スチューデントのt分布(Student's t-Distribution)のグラフを作成します。t分布については「1次元スチューデントのt分布の定義式の確認 - からっぽのしょこ」を参照してください。
利用するパッケージを読み込みます。
# 利用パッケージ library(tidyverse) library(gganimate) library(LaplacesDemon)
この記事では、基本的にパッケージ名::関数名()の記法を使うので、パッケージを読み込む必要はありません。ただし、作図コードがごちゃごちゃしないようにパッケージ名を省略しているため、ggplot2は読み込む必要があります。
magrittrパッケージのパイプ演算子%>%ではなく、ベースパイプ(ネイティブパイプ)演算子|>を使っています。%>%に置き換えても処理できます。
分布の変化をアニメーション(gif画像)で確認するのにgganimateパッケージを利用します。不要であれば省略してください。
一般化したt分布の計算にLaplacesDemonパッケージを利用します。
定義式の確認
まずは、標準化されたスチューデントのt分布と一般化された(標準化されていない)スチューデントのt分布の定義式を確認します。(標準化・一般化という表現は雰囲気です。正しい言い方があれば教えてください。)
標準化t分布は、次の式で定義されます。
ここで、$\nu$は自由度、$\pi$は円周率、$\Gamma(x)$はガンマ関数です。自由度は形状パラメータとも呼ばれます。
確率変数の実現値$x$は、実数となります。$\nu$は、正の整数を満たす必要があります。
また、一般化t分布は、次の2つの式で定義されます。
ここで、$\mu$は位置パラメータ、$\sigma$はスケールパラメータ、$\lambda$は逆スケールパラメータです。
$\mu$は実数、$\sigma$は正の実数、$\lambda$は正の実数を満たす必要があります。また、$\sigma = \frac{1}{\sqrt{\lambda}}$、$\lambda = \frac{1}{\sigma^2}$の関係が成り立ちます。
$\mu = 0$、$\sigma = \lambda = 1$のとき、標準化t分布と一般化t分布は一致します。
t分布の期待値・分散・最頻値は、次の式で計算できます。詳しくはいつか書きます。
これらの計算を行いグラフを作成します。
グラフの作成
ggplot2パッケージを利用して、t分布のグラフを作成します。t分布の確率や統計量の計算については「【R】1次元スチューデントのt分布の計算 - からっぽのしょこ」を参照してください。
t分布のパラメータ$\nu, \mu, \sigma$を設定します。
# 形状パラメータ(自由度)を指定 nu <- 5 # 位置パラメータを指定 mu <- 2 # スケールパラメータを指定 sigma <- 0.5
形状パラメータ(自由度)(正の整数)$\nu$、位置パラメータ(実数)$\mu$、スケールパラメータ(正の実数)$\sigma$を指定します。
t分布の確率変数がとり得る値$x$を作成します。
# xの値を作成 x_vals <- seq(from = mu - sigma*5, to = mu + sigma*5, length.out = 251)
$x$の値を作成してx_valsとします。この例では、muを中心にスケールパラメータの5倍を範囲とします。
$x$の値ごとに確率密度を計算します。
# スチューデントのt分布を計算 dens_df <- tibble::tibble( x = x_vals, # 確率変数 density = LaplacesDemon::dst(x = x_vals, mu = mu, sigma = sigma, nu = nu) # 確率密度 ) dens_df
## # A tibble: 251 × 2 ## x density ## <dbl> <dbl> ## 1 -0.5 0.00351 ## 2 -0.48 0.00366 ## 3 -0.46 0.00381 ## 4 -0.44 0.00397 ## 5 -0.42 0.00413 ## 6 -0.4 0.00430 ## 7 -0.38 0.00449 ## 8 -0.36 0.00468 ## 9 -0.34 0.00487 ## 10 -0.32 0.00508 ## # … with 241 more rows
t分布の確率密度は、LaplacesDemonパッケージのdst()で計算できます。確率変数の引数x
にx_vals、位置パラメータの引数muにmu、スケールパラメータの引数sigmaにsigma、形状パラメータ(自由度)の引数nuにnuを指定します。逆スケールパラメータ$\lambda$を使う場合は、sigma引数に$\lambda$の逆数の平方根$\frac{1}{\sqrt{\lambda}}$を指定します。
t分布のグラフを作成します。
# スチューデントのt分布のグラフを作成 ggplot(data = dens_df, mapping = aes(x = x, y = density)) + # データ geom_line(color = "#00A968", size = 1) + # 折れ線グラフ labs( title = "Student's t Distribution", subtitle = paste0("nu=", nu, ", mu=", mu, ", sigma=", sigma), # (文字列表記用) #subtitle = parse(text = paste0("list(nu==", nu, ", mu==", mu, ", sigma==", sigma, ")")), # (数式表記用) x = "x", y = "density" ) # ラベル
ギリシャ文字などの記号を使った数式を表示する場合は、expression()の記法を使います。等号は"=="、複数の(数式上の)変数を並べるには"list(変数1, 変数2)"とします。(プログラム上の)変数の値を使う場合は、parse()のtext引数に指定します。
この分布に統計量の情報を重ねて表示します。
# 補助線用の統計量を計算 E_x <- mu # (nu > 1) s_x <- sqrt(sigma^2 * nu / (nu - 2.0)) # (nu > 2) # 統計量を重ねたt分布のグラフを作成:線のみ ggplot(data = dens_df, mapping = aes(x = x, y = density)) + # データ geom_line(color = "#00A968", size = 1) + # 分布 geom_vline(xintercept = E_x, color = "blue", size = 1, linetype = "dashed") + # 期待値 geom_vline(xintercept = E_x-s_x, color = "orange", size = 1, linetype = "dotted") + # 期待値 - 標準偏差 geom_vline(xintercept = E_x+s_x, color = "orange", size = 1, linetype = "dotted") + # 期待値 + 標準偏差 labs( title = "Student's t Distribution", subtitle = parse(text = paste0("list(nu==", nu, ", mu==", mu, ", sigma==", sigma, ")")), x = "x", y = "density" ) # ラベル
期待値・標準偏差(分散の平方根)を計算して、それぞれgeom_vline()で垂直線を引きます。標準偏差については期待値から足し引きした値を使います。
凡例を表示する場合は、垂直線を引く値をデータフレームにまとめます。
# 統計量を格納 stat_df <- tibble::tibble( statistic = c(mu, mu-sigma, mu+sigma), # 統計量 type = c("mean", "sd", "sd") # 色分け用ラベル ) stat_df
## # A tibble: 3 × 2 ## statistic type ## <dbl> <chr> ## 1 2 mean ## 2 1.5 sd ## 3 2.5 sd
各統計量を区別するための文字列も格納します。期待値と標準偏差の差・和は同じラベルとします。
線の色や種類、凡例テキストを指定する場合は、設定用の名前付きベクトルを作成します。
# 凡例用の設定を作成:(数式表示用) color_vec <- c(mean = "blue", sd = "orange") linetype_vec <- c(mean = "dashed", sd = "dotted") label_vec <- c(mean = expression(E(x)), sd = expression(E(x) %+-% sqrt(V(x)))) color_vec; linetype_vec; label_vec
## mean sd ## "blue" "orange" ## mean sd ## "dashed" "dotted" ## expression(mean = E(x), sd = E(x) %+-% sqrt(V(x)))
各要素(設定)の名前をtype列の文字列と対応させて、設定(線の色・種類と数式用のテキスト)を指定します。
凡例付きのグラフを作成します。
# 統計量を重ねたt分布のグラフを作成:凡例付き ggplot() + geom_line(data = dens_df, mapping = aes(x = x, y = density), color = "#00A968", size = 1) + # 分布 geom_vline(data = stat_df, mapping = aes(xintercept = statistic, color = type, linetype = type), size = 1) + # 期待値 scale_color_manual(values = color_vec, labels = label_vec, name = "statistic") + # 線の色:(色指定と数式表示用) scale_linetype_manual(values = linetype_vec, labels = label_vec, name = "statistic") + # 線の種類:(線指定と数式表示用) guides(color = guide_legend(override.aes = list(size = 0.5))) + # 凡例の体裁 theme(legend.text.align = 0) + # 図の体裁:凡例 labs(title = "Student's t Distribution", subtitle = parse(text = paste0("list(nu==", nu, ", mu==", mu, ", sigma==", sigma, ")")), x = "x", y = "density") # ラベル
scale_color_manual()のvalues引数に線の色、scale_linetype_manual()のvalues引数に線の種類を指定します。凡例テキストを指定する場合は、それぞれのlabels引数に指定します。names引数は凡例ラベルです。
凡例テキストは、theme()のlegend.text.align引数に0を指定すると左寄せ、デフォルト(1)だと右寄せになります。
$\nu > 1$のとき期待値と最頻値が一致します。$\nu \leq 0$のときは期待値を計算できません。
t分布のグラフを描画できました。以降は、ここまでの作図処理を用いて、パラメータの影響を確認していきます。
パラメータと分布の関係を並べて比較
複数のパラメータのグラフを比較することで、パラメータの値と分布の形状の関係を確認します。
自由度の影響
複数の$\nu$を指定し、$\mu, \sigma$を固定して、それぞれ分布を作図します。
・作図コード(クリックで展開)
# 自由度として利用する値を指定 nu_vals <- c(1, 2, 3, 4, 5, 10, 100) # 固定するパラメータを指定 mu <- 2 sigma <- 0.5 # xの値を作成 x_vals <- seq(from = mu - sigma*5, to = mu + sigma*5, length.out = 251) # パラメータごとにt分布を計算 res_dens_df <- tidyr::expand_grid( x = x_vals, nu = nu_vals ) |> # 全ての組み合わせを作成 dplyr::arrange(nu, x) |> # パラメータごとに並べ替え dplyr::mutate( density = LaplacesDemon::dst(x = x_vals, mu = mu, sigma = sigma, nu = nu), parameter = paste0("nu=", nu, ", mu=", mu, ", sigma=", sigma) |> factor(levels = paste0("nu=", sort(nu_vals), ", mu=", mu, ", sigma=", sigma)) # 色分け用ラベル ) # 確率密度を計算 res_dens_df
## # A tibble: 1,757 × 4 ## x nu density parameter ## <dbl> <dbl> <dbl> <fct> ## 1 -0.5 1 0.0245 nu=1, mu=2, sigma=0.5 ## 2 -0.48 1 0.0249 nu=1, mu=2, sigma=0.5 ## 3 -0.46 1 0.0253 nu=1, mu=2, sigma=0.5 ## 4 -0.44 1 0.0257 nu=1, mu=2, sigma=0.5 ## 5 -0.42 1 0.0261 nu=1, mu=2, sigma=0.5 ## 6 -0.4 1 0.0265 nu=1, mu=2, sigma=0.5 ## 7 -0.38 1 0.0269 nu=1, mu=2, sigma=0.5 ## 8 -0.36 1 0.0273 nu=1, mu=2, sigma=0.5 ## 9 -0.34 1 0.0278 nu=1, mu=2, sigma=0.5 ## 10 -0.32 1 0.0283 nu=1, mu=2, sigma=0.5 ## # … with 1,747 more rows
確率変数x_valsとパラメータnu_valsそれぞれの要素の全ての組み合わせをexpand_grid()で作成します。
組み合わせごとに確率密度を計算して、パラメータごとにラベルを作成します。ラベルは、グラフの設定や凡例テキストとして使います。文字列型だと文字列の基準で順序が決まるので、因子型にしてパラメータに応じたレベル(順序)を設定します。
凡例を数式で表示する場合は、expression()に対応した記法に変換します。
# 凡例用のラベルを作成:(数式表示用) label_vec <- res_dens_df[["parameter"]] |> stringr::str_replace_all(pattern = "=", replacement = "==") |> # 等号表示用の記法に変換 (\(.){paste0("list(", ., ")")})() |> # カンマ表示用の記法に変換 unique() |> # 重複を除去 parse(text = _) # expression関数化 names(label_vec) <- unique(res_dens_df[["parameter"]]) # ggplotに指定する文字列に対応する名前付きベクトルに変換 label_vec[1]
## expression(`nu=1, mu=2, sigma=0.5` = list(nu == 1, mu == 2, sigma == ## 0.5))
変換後の文字列のベクトルに対してnames()を使って、元の文字列を各要素の名前として設定します。
3行目の処理では、無名関数function()の省略記法\()を使って、(\(引数){引数を使った具体的な処理})()としています。直前のパイプ演算子を%>%にすると、行全体(\(引数){処理})()を{}の中の処理(この例だとpaste0("list(", ., ")"))に置き換えられます(置き換えられるように引数名を.にしています)。
res_prob_dfのparameter列の要素nu=1, mu=2, sigma=0.5と、変換後の文字列list(nu == 1, mu == 2, sigma == 0.5)が対応しています。
パラメータごとにグラフを作成します。
# パラメータごとにt分布を作図 ggplot(data = res_dens_df, mapping = aes(x = x, y = density, color = parameter)) + # データ geom_line(size = 1) + # 折れ線グラフ scale_color_hue(labels = label_vec) + # 線の色:(数式表示用) theme(legend.text.align = 0) + # 図の体裁:凡例 labs(title = "Student's t Distribution", x = "x", y = "density") # ラベル
$\nu$が小さいほど裾が厚くなり、$\nu$が大きいほどガウス分布に近付きます。
位置パラメータの影響
複数の$\mu$を指定し、$\nu, \sigma$を固定して、それぞれ分布を作図します。
・作図コード(クリックで展開)
# 位置パラメータとして利用する値を指定 mu_vals <- c(-3.5, -2, -0.5, 1, 2.5) # 固定するパラメータを指定 nu <- 1 sigma <- 0.5 # xの値を作成 x_vals <- seq(from = min(mu_vals)-sigma, to = max(mu_vals)+sigma, length.out = 251) # パラメータごとにt分布を計算 res_dens_df <- tidyr::expand_grid( x = x_vals, mu = mu_vals ) |> # 全ての組み合わせを作成 dplyr::arrange(mu, x) |> # パラメータごとに並べ替え dplyr::mutate( density = LaplacesDemon::dst(x = x_vals, mu = mu, sigma = sigma, nu = nu), parameter = paste0("nu=", nu, ", mu=", mu, ", sigma=", sigma) |> factor(levels = paste0("nu=", nu, ", mu=", sort(mu_vals), ", sigma=", sigma)) # 色分け用ラベル ) # 確率密度を計算 res_dens_df
## # A tibble: 1,255 × 4 ## x mu density parameter ## <dbl> <dbl> <dbl> <fct> ## 1 -4 -3.5 0.318 nu=1, mu=-3.5, sigma=0.5 ## 2 -3.97 -3.5 0.337 nu=1, mu=-3.5, sigma=0.5 ## 3 -3.94 -3.5 0.356 nu=1, mu=-3.5, sigma=0.5 ## 4 -3.92 -3.5 0.376 nu=1, mu=-3.5, sigma=0.5 ## 5 -3.89 -3.5 0.397 nu=1, mu=-3.5, sigma=0.5 ## 6 -3.86 -3.5 0.419 nu=1, mu=-3.5, sigma=0.5 ## 7 -3.83 -3.5 0.442 nu=1, mu=-3.5, sigma=0.5 ## 8 -3.80 -3.5 0.465 nu=1, mu=-3.5, sigma=0.5 ## 9 -3.78 -3.5 0.488 nu=1, mu=-3.5, sigma=0.5 ## 10 -3.75 -3.5 0.511 nu=1, mu=-3.5, sigma=0.5 ## # … with 1,245 more rows
x_valsとmu_valsの全ての組み合わせを作成して、同様に処理します。
作図については同じコードで処理できます。
$\mu$が大きいほど、山が右に移動します。
スケールパラメータの影響
複数の$\sigma$を指定し、$\nu, \mu$を固定して、それぞれ分布を作図します。
・作図コード(クリックで展開)
# スケールパラメータとして利用する値を指定 sigma_vals <- c(0.5, 1, 2, 4, 8) # 固定するパラメータを指定 nu <- 1 mu <- 2 # xの値を作成 x_vals <- seq(from = mu - max(sigma_vals), to = mu + max(sigma_vals), length.out = 251) # パラメータごとにt分布を計算 res_dens_df <- tidyr::expand_grid( x = x_vals, sigma = sigma_vals ) |> # 全ての組み合わせを作成 dplyr::arrange(sigma, x) |> # パラメータごとに並べ替え dplyr::mutate( density = LaplacesDemon::dst(x = x_vals, mu = mu, sigma = sigma, nu = nu), parameter = paste0("nu=", nu, ", mu=", mu, ", sigma=", sigma) |> factor(levels = paste0("nu=", nu, ", mu=", mu, ", sigma=", sort(sigma_vals))) # 色分け用ラベル ) # 確率密度を計算 res_dens_df
## # A tibble: 1,255 × 4 ## x sigma density parameter ## <dbl> <dbl> <dbl> <fct> ## 1 -6 0.5 0.00248 nu=1, mu=2, sigma=0.5 ## 2 -5.94 0.5 0.00252 nu=1, mu=2, sigma=0.5 ## 3 -5.87 0.5 0.00256 nu=1, mu=2, sigma=0.5 ## 4 -5.81 0.5 0.00260 nu=1, mu=2, sigma=0.5 ## 5 -5.74 0.5 0.00264 nu=1, mu=2, sigma=0.5 ## 6 -5.68 0.5 0.00269 nu=1, mu=2, sigma=0.5 ## 7 -5.62 0.5 0.00273 nu=1, mu=2, sigma=0.5 ## 8 -5.55 0.5 0.00278 nu=1, mu=2, sigma=0.5 ## 9 -5.49 0.5 0.00283 nu=1, mu=2, sigma=0.5 ## 10 -5.42 0.5 0.00287 nu=1, mu=2, sigma=0.5 ## # … with 1,245 more rows
x_valsとsigma_valsの全ての組み合わせを作成して、同様に処理します。
作図については同じコードで処理できます。
$\sigma$が小さいほど、山が細く高くなります。
逆スケールパラメータの影響
$\sigma$の代わりに複数の$\lambda$を指定し、$\nu, \mu$を固定して、それぞれ分布を作図します。
・作図コード(クリックで展開)
# 逆スケールパラメータとして利用する値を指定 lambda_vals <- c(0.25, 0.5, 1, 2, 4) # 固定するパラメータを指定 nu <- 1 mu <- 2 # xの値を作成 sigma <- 1 / sqrt(min(lambda_vals)) x_vals <- seq(from = mu - sigma*4, to = mu + sigma*4, length.out = 251) # パラメータごとにt分布を計算 res_dens_df <- tidyr::expand_grid( x = x_vals, lambda = lambda_vals ) |> # 全ての組み合わせを作成 dplyr::arrange(lambda, x) |> # パラメータごとに並べ替え dplyr::mutate( density = LaplacesDemon::dst(x = x_vals, mu = mu, sigma = 1/sqrt(lambda), nu = nu), parameter = paste0("nu=", nu, ", mu=", mu, ", lambda=", lambda) |> factor(levels = paste0("nu=", nu, ", mu=", mu, ", lambda=", sort(lambda_vals))) # 色分け用ラベル ) # 確率密度を計算 res_dens_df
## # A tibble: 1,255 × 4 ## x lambda density parameter ## <dbl> <dbl> <dbl> <fct> ## 1 -6 0.25 0.00936 nu=1, mu=2, lambda=0.25 ## 2 -5.94 0.25 0.00950 nu=1, mu=2, lambda=0.25 ## 3 -5.87 0.25 0.00965 nu=1, mu=2, lambda=0.25 ## 4 -5.81 0.25 0.00980 nu=1, mu=2, lambda=0.25 ## 5 -5.74 0.25 0.00995 nu=1, mu=2, lambda=0.25 ## 6 -5.68 0.25 0.0101 nu=1, mu=2, lambda=0.25 ## 7 -5.62 0.25 0.0103 nu=1, mu=2, lambda=0.25 ## 8 -5.55 0.25 0.0104 nu=1, mu=2, lambda=0.25 ## 9 -5.49 0.25 0.0106 nu=1, mu=2, lambda=0.25 ## 10 -5.42 0.25 0.0108 nu=1, mu=2, lambda=0.25 ## # … with 1,245 more rows
x_valsとlambda_valsの全ての組み合わせを作成して、同様に処理します。
dst()のsigma引数にlambda_valsの値の逆数の平方根を指定します。
作図については同じコードで処理できます。
$\lambda$が大きいほど、山が細く高くなります。
パラメータと分布の関係をアニメーションで可視化
前節では、複数のパラメータのグラフを並べて比較しました。次は、パラメータの値を少しずつ変化させて、分布の形状の変化をアニメーションで確認します。
自由度の影響
$\nu$の値を変化させ、$\mu, \sigma$を固定して、それぞれ分布を作図します。
・作図コード(クリックで展開)
# 自由度として利用する値を指定 nu_vals <- 1:50 # (線の変化用) #nu_vals <- seq(from = 1, to = 25, by = 1) # (線の軌跡用) length(nu_vals) # フレーム数 # 固定するパラメータを指定 mu <- 2 sigma <- 0.5 # xの値を作成 x_vals <- seq(from = mu - sigma*5, to = mu + sigma*5, length.out = 251) # パラメータごとにt分布を計算 anime_dens_df <- tidyr::expand_grid( x = x_vals, nu = nu_vals ) |> # 全ての組み合わせを作成 dplyr::arrange(nu, x) |> # パラメータごとに並べ替え dplyr::mutate( density = LaplacesDemon::dst(x = x_vals, mu = mu, sigma = sigma, nu = nu), parameter = paste0("nu=", round(nu, 2), ", mu=", mu, ", sigma=", sigma) |> factor(levels = paste0("nu=", round(nu_vals, 2), ", mu=", mu, ", sigma=", sigma)) # 色分け用ラベル ) # 確率密度を計算 anime_dens_df
## # A tibble: 12,550 × 4 ## x nu density parameter ## <dbl> <int> <dbl> <fct> ## 1 -0.5 1 0.0245 nu=1, mu=2, sigma=0.5 ## 2 -0.48 1 0.0249 nu=1, mu=2, sigma=0.5 ## 3 -0.46 1 0.0253 nu=1, mu=2, sigma=0.5 ## 4 -0.44 1 0.0257 nu=1, mu=2, sigma=0.5 ## 5 -0.42 1 0.0261 nu=1, mu=2, sigma=0.5 ## 6 -0.4 1 0.0265 nu=1, mu=2, sigma=0.5 ## 7 -0.38 1 0.0269 nu=1, mu=2, sigma=0.5 ## 8 -0.36 1 0.0273 nu=1, mu=2, sigma=0.5 ## 9 -0.34 1 0.0278 nu=1, mu=2, sigma=0.5 ## 10 -0.32 1 0.0283 nu=1, mu=2, sigma=0.5 ## # … with 12,540 more rows
値の間隔が一定になるようにnu_valsを作成します。パラメータごとにフレームを切り替えるので、mu_valsの要素数がアニメーションのフレーム数になります。
データフレームの変数名以外は「並べて比較」のときと同じ処理です。
スチューデントのt分布のアニメーション(gif画像)を作成します。
# t分布のアニメーションを作図 anime_dens_graph <- ggplot(data = anime_dens_df, mapping = aes(x = x, y = density)) + # データ geom_line(color = "#00A968", size = 1) + # 折れ線グラフ gganimate::transition_manual(parameter) + # フレーム labs(title = "Student's t Distribution", subtitle = "{current_frame}", x = "x", y = "density") # ラベル # gif画像を作成 gganimate::animate(anime_dens_graph, nframes = length(nu_vals), fps = 10, width = 800, height = 600)
transition_manual()にフレームの順序を表す列を指定します。この例では、因子型のラベルのレベルの順に描画されます。
animate()のフレーム数の引数nframesにパラメータ数、フレームレートの引数fpsに1秒当たりのフレーム数を指定します。fps引数の値が大きいほどフレームが早く切り替わります。ただし、値が大きいと意図通りに動作しません。
続いて、分布の変化の軌跡を描画するアニメーションを作成します。mu_valsの値の間隔を粗くしておきます。
# t分布のアニメーションを作図 anime_dens_graph <- ggplot(data = anime_dens_df, mapping = aes(x = x, y = density, color = factor(nu), group = nu)) + # データ geom_line() + # 折れ線グラフ gganimate::transition_reveal(nu) + # フレーム labs(title = "Student's t Distribution", subtitle = paste0("nu={frame_along}, mu=", mu, ", sigma=", sigma), color = expression(nu), x = "x", y = "density") # ラベル # gif画像を作成 gganimate::animate(anime_dens_graph, nframes = length(nu_vals)+10, end_pause = 10, fps = 10, width = 800, height = 600)
transition_reveal()を使ってフレームを切り替えると、過去フレームのグラフが表示され続けます。
animate()のend_pause引数を指定すると、最後のフレームでグラフが一時停止(最後のグラフを指定したフレーム数表示)します。
位置パラメータの影響
続いて、$\mu$の値を変化させ、$\nu, \sigma$を固定して、それぞれ分布を計算します。
・作図コード(クリックで展開)
# 位置パラメータとして利用する値を指定 mu_vals <- seq(from = -3, to = 3, by = 0.1) # (線の変化用) #mu_vals <- seq(from = -3, to = 3, by = 0.5) # (線の軌跡用) length(mu_vals) # フレーム数 # 固定するパラメータを指定 nu <- 1 sigma <- 0.5 # xの値を作成 x_vals <- seq(from = min(mu_vals)-sigma, to = max(mu_vals)+sigma, length.out = 251) # パラメータごとにt分布を計算 anime_dens_df <- tidyr::expand_grid( x = x_vals, mu = mu_vals ) |> # 全ての組み合わせを作成 dplyr::arrange(mu, x) |> # パラメータごとに並べ替え dplyr::mutate( density = LaplacesDemon::dst(x = x_vals, mu = mu, sigma = sigma, nu = nu), parameter = paste0("nu=", nu, ", mu=", round(mu, 2), ", sigma=", sigma) |> factor(levels = paste0("nu=", nu, ", mu=", round(mu_vals, 2), ", sigma=", sigma)) # 色分け用ラベル ) # 確率密度を計算 anime_dens_df
## # A tibble: 15,311 × 4 ## x mu density parameter ## <dbl> <dbl> <dbl> <fct> ## 1 -3.5 -3 0.318 nu=1, mu=-3, sigma=0.5 ## 2 -3.47 -3 0.337 nu=1, mu=-3, sigma=0.5 ## 3 -3.44 -3 0.356 nu=1, mu=-3, sigma=0.5 ## 4 -3.42 -3 0.376 nu=1, mu=-3, sigma=0.5 ## 5 -3.39 -3 0.397 nu=1, mu=-3, sigma=0.5 ## 6 -3.36 -3 0.419 nu=1, mu=-3, sigma=0.5 ## 7 -3.33 -3 0.442 nu=1, mu=-3, sigma=0.5 ## 8 -3.30 -3 0.465 nu=1, mu=-3, sigma=0.5 ## 9 -3.28 -3 0.488 nu=1, mu=-3, sigma=0.5 ## 10 -3.25 -3 0.511 nu=1, mu=-3, sigma=0.5 ## # … with 15,301 more rows
「並べて比較」や「自由度の影響」と同様に処理します。
「自由度の影響」のコードで作図できます。フレーム数はmu_valsの要素数です。
# gif画像を作成 gganimate::animate(anime_dens_graph, nframes = length(mu_vals), fps = 10, width = 800, height = 600)
続いて、分布の変化の軌跡を描画するアニメーションを作成します。
# t分布のアニメーションを作図 anime_dens_graph <- ggplot(data = anime_dens_df, mapping = aes(x = x, y = density, color = factor(mu), group = mu)) + # データ geom_line() + # 折れ線グラフ gganimate::transition_reveal(mu) + # フレーム labs(title = "Student's t Distribution", subtitle = paste0("nu=", nu, ", mu={frame_along}, sigma=", sigma), color = expression(mu), x = "x", y = "density") # ラベル # gif画像を作成 gganimate::animate(anime_dens_graph, nframes = length(mu_vals)+10, end_pause = 10, fps = 10, width = 800, height = 600)
「自由度の影響」のnuをmuに置き換えて処理します。
スケールパラメータの影響
$\sigma$の値を変化させ、$\nu, \mu$を固定して、それぞれ分布を計算します。
・作図コード(クリックで展開)
# スケールパラメータとして利用する値を指定 sigma_vals <- seq(from = 0.5, to = 5, by = 0.1) # (線の変化用) #sigma_vals <- seq(from = 0.5, to = 10, by = 0.5) # (線の軌跡用) length(sigma_vals) # フレーム数 # 固定するパラメータを指定 nu <- 1 mu <- 2 # xの値を作成 x_vals <- seq(from = mu - max(sigma_vals), to = mu + max(sigma_vals), length.out = 251) # パラメータごとにt分布を計算 anime_dens_df <- tidyr::expand_grid( x = x_vals, sigma = sigma_vals ) |> # 全ての組み合わせを作成 dplyr::arrange(sigma, x) |> # パラメータごとに並べ替え dplyr::mutate( density = LaplacesDemon::dst(x = x_vals, mu = mu, sigma = sigma, nu = nu), parameter = paste0("nu=", nu, ", mu=", mu, ", sigma=", round(sigma, 2)) |> factor(levels = paste0("nu=", nu, ", mu=", mu, ", sigma=", round(sigma_vals, 2))) # 色分け用ラベル ) # 確率密度を計算 anime_dens_df
## # A tibble: 11,546 × 4 ## x sigma density parameter ## <dbl> <dbl> <dbl> <fct> ## 1 -3 0.5 0.00630 nu=1, mu=2, sigma=0.5 ## 2 -2.96 0.5 0.00640 nu=1, mu=2, sigma=0.5 ## 3 -2.92 0.5 0.00651 nu=1, mu=2, sigma=0.5 ## 4 -2.88 0.5 0.00661 nu=1, mu=2, sigma=0.5 ## 5 -2.84 0.5 0.00672 nu=1, mu=2, sigma=0.5 ## 6 -2.8 0.5 0.00683 nu=1, mu=2, sigma=0.5 ## 7 -2.76 0.5 0.00695 nu=1, mu=2, sigma=0.5 ## 8 -2.72 0.5 0.00706 nu=1, mu=2, sigma=0.5 ## 9 -2.68 0.5 0.00718 nu=1, mu=2, sigma=0.5 ## 10 -2.64 0.5 0.00731 nu=1, mu=2, sigma=0.5 ## # … with 11,536 more rows
「並べて比較」や「自由度の影響」と同様に処理します。
「自由度の影響」のコードで作図できます。フレーム数はsigma_valsの要素数です。
# gif画像を作成 gganimate::animate(anime_dens_graph, nframes = length(sigma_vals), fps = 10, width = 800, height = 600)
続いて、分布の変化の軌跡を描画するアニメーションを作成します。
# t分布のアニメーションを作図 anime_dens_graph <- ggplot(data = anime_dens_df, mapping = aes(x = x, y = density, color = factor(sigma), group = sigma)) + # データ geom_line() + # 折れ線グラフ gganimate::transition_reveal(sigma) + # フレーム labs(title = "Student's t Distribution", subtitle = paste0("nu=", nu, ", mu=", mu, ", sigma={frame_along}"), color = expression(sigma), x = "x", y = "density") # ラベル # gif画像を作成 gganimate::animate(anime_dens_graph, nframes = length(sigma_vals)+10, end_pause=10, fps = 10, width = 800, height = 600)
「自由度の影響」のnuをsigmaに置き換えて処理します。
逆スケールパラメータの影響
$\sigma$の代わりに$\lambda$の値を変化させ、$\nu, \mu$を固定して、それぞれ分布を計算します。
・作図コード(クリックで展開)
# 逆スケールパラメータとして利用する値を指定 lambda_vals <- seq(from = 0.1, to = 10, by = 0.1) # (線の変化用) #lambda_vals <- seq(from = 0.5, to = 10, by = 0.5) # (線の軌跡用) length(lambda_vals) # フレーム数 # 固定するパラメータを指定 nu <- 1 mu <- 2 # xの値を作成 sigma <- 1 / sqrt(min(lambda_vals)) x_vals <- seq(from = mu - sigma*2, to = mu + sigma*2, length.out = 251) # パラメータごとにt分布を計算 anime_dens_df <- tidyr::expand_grid( x = x_vals, lambda = lambda_vals ) |> # 全ての組み合わせを作成 dplyr::arrange(lambda, x) |> # パラメータごとに並べ替え dplyr::mutate( density = LaplacesDemon::dst(x = x_vals, mu = mu, sigma = 1/sqrt(lambda), nu = nu), parameter = paste0("nu=", nu, ", mu=", mu, ", lambda=", round(lambda, 2)) |> factor(levels = paste0("nu=", nu, ", mu=", mu, ", lambda=", round(lambda_vals, 2))) # 色分け用ラベル ) # 確率密度を計算 anime_dens_df
## # A tibble: 25,100 × 4 ## x lambda density parameter ## <dbl> <dbl> <dbl> <fct> ## 1 -4.32 0.1 0.0201 nu=1, mu=2, lambda=0.1 ## 2 -4.27 0.1 0.0204 nu=1, mu=2, lambda=0.1 ## 3 -4.22 0.1 0.0207 nu=1, mu=2, lambda=0.1 ## 4 -4.17 0.1 0.0209 nu=1, mu=2, lambda=0.1 ## 5 -4.12 0.1 0.0212 nu=1, mu=2, lambda=0.1 ## 6 -4.07 0.1 0.0215 nu=1, mu=2, lambda=0.1 ## 7 -4.02 0.1 0.0218 nu=1, mu=2, lambda=0.1 ## 8 -3.97 0.1 0.0221 nu=1, mu=2, lambda=0.1 ## 9 -3.92 0.1 0.0223 nu=1, mu=2, lambda=0.1 ## 10 -3.87 0.1 0.0226 nu=1, mu=2, lambda=0.1 ## # … with 25,090 more rows
「並べて比較」や「自由度の影響」と同様に処理します。
「自由度の影響」のコードで作図できます。フレーム数はlambda_valsの要素数です。
# gif画像を作成 gganimate::animate(anime_dens_graph, nframes = length(lambda_vals), fps = 10, width = 800, height = 600)
続いて、分布の変化の軌跡を描画するアニメーションを作成します。
# t分布のアニメーションを作図 anime_dens_graph <- ggplot(data = anime_dens_df, mapping = aes(x = x, y = density, color = factor(lambda), group = lambda)) + # データ geom_line() + # 折れ線グラフ gganimate::transition_reveal(lambda) + # フレーム labs(title = "Student's t Distribution", subtitle = paste0("nu=", nu, ", mu=", mu, ", lambda={frame_along}"), color = expression(lambda), x = "x", y = "density") # ラベル # gif画像を作成 gganimate::animate(anime_dens_graph, nframes = length(lambda_vals)+10, end_pause=10, fps = 10, width = 800, height = 600)
「自由度の影響」のnuをlambdaに置き換えて処理します。
パラメータと統計量の関係をアニメーションで可視化
ここまでは、パラメータと分布の関係を確認しました。次は、パラメータと統計量の関係をアニメーションで確認します。
自由度の影響
「パラメータと分布の形状の関係」と同様に、$\nu$の値を変化させ、$\mu, \sigma$を固定して、anime_dens_dfを作成します。
・作図コード(クリックで展開)
パラメータごとに統計量を計算します。
# パラメータごとに統計量を計算 anime_stat_df <- tibble::tibble( nu = nu_vals, mean = dplyr::if_else( nu > 1, true = mu, false = as.numeric(NA) ), # 期待値 sd = dplyr::if_else( nu > 2, true = sqrt(sigma^2 * nu / (nu - 2)), false = as.numeric(NA) ), # 標準偏差 parameter = paste0("nu=", round(nu_vals, 2), ", mu=", mu, ", sigma=", sigma) |> factor(levels = paste0("nu=", round(nu_vals, 2), ", mu=", mu, ", sigma=", sigma)) # フレーム切替用のラベル ) |> # 統計量を計算 dplyr::mutate( sd_minus = mean - sd, sd_plus = mean + sd ) |> # 期待値±標準偏差を計算 dplyr::select(!c(nu, sd)) |> # 不要な列を削除 tidyr::pivot_longer( cols = !parameter, names_to = "type", values_to = "statistic" ) |> # 統計量の列をまとめる dplyr::mutate( type = stringr::str_replace(type, pattern = "sd_.*", replacement = "sd")) # 期待値±標準偏差のカテゴリを統一 anime_stat_df
## # A tibble: 150 × 3 ## parameter type statistic ## <fct> <chr> <dbl> ## 1 nu=1, mu=2, sigma=0.5 mean NA ## 2 nu=1, mu=2, sigma=0.5 sd NA ## 3 nu=1, mu=2, sigma=0.5 sd NA ## 4 nu=2, mu=2, sigma=0.5 mean 2 ## 5 nu=2, mu=2, sigma=0.5 sd NA ## 6 nu=2, mu=2, sigma=0.5 sd NA ## 7 nu=3, mu=2, sigma=0.5 mean 2 ## 8 nu=3, mu=2, sigma=0.5 sd 1.13 ## 9 nu=3, mu=2, sigma=0.5 sd 2.87 ## 10 nu=4, mu=2, sigma=0.5 mean 2 ## # … with 140 more rows
期待値・標準偏差・最頻値を計算して、さらに期待値から標準偏差を引いた値と足した値を求めます。
期待値・標準偏差の和と差・最頻値の4つの列をpivot_longer()でまとめます。列名がラベルになるので、標準偏差の和と差のラベルを統一します。
線と凡例の設定用の名前付きベクトルを作成します。
# 凡例用の設定を作成:(数式表示用) color_vec <- c(mean = "blue", sd = "orange") linetype_vec <- c(mean = "dashed", sd = "dotted") label_vec <- c(mean = expression(E(x)), sd = expression(E(x) %+-% sqrt(V(x))))
「グラフの作成」のときと同じ処理です。
統計量の情報を重ねたガンマ分布のアニメーション(gif画像)を作成します。
# 統計量を重ねたt分布のアニメーションを作図 anime_dens_graph <- ggplot() + geom_line(data = anime_dens_df, mapping = aes(x = x, y = density), color = "#00A968", size = 1) + # 分布 geom_vline(data = anime_stat_df, mapping = aes(xintercept = statistic, color = type, linetype = type), size = 1) + # 統計量 gganimate::transition_manual(parameter) + # フレーム scale_color_manual(values = color_vec, labels = label_vec, name = "statistic") + # 線の色:(色指定と数式表示用) scale_linetype_manual(values = linetype_vec, labels = label_vec, name = "statistic") + # 線の種類:(線指定と数式表示用) guides(color = guide_legend(override.aes = list(size = 0.5))) + # 凡例の体裁 theme(legend.text.align = 0) + # 図の体裁:凡例 #coord_cartesian(xlim = c(min(x_vals), max(x_vals))) + # 軸の表示範囲 labs(title = "Student's t Distribution", subtitle = "{current_frame}", x = "x", y = "density") # ラベル # gif画像を作成 gganimate::animate(anime_dens_graph, nframes = length(nu_vals), fps = 10, width = 800, height = 600)
位置パラメータの影響
$\mu$の値を変化させ、$\nu, \sigma$を固定して、anime_dens_dfを作成します。
・作図コード(クリックで展開)
パラメータごとに統計量を計算します。
# パラメータごとに統計量を計算 anime_stat_df <- tibble::tibble( nu = nu, mu = mu_vals, mean = dplyr::if_else( nu > 1, true = mu, false = as.numeric(NA) ), # 期待値 sd = dplyr::if_else( nu > 2, true = sqrt(sigma^2 * nu / (nu - 2)), false = as.numeric(NA) ), # 標準偏差 parameter = paste0("nu=", nu, ", mu=", round(mu_vals, 2), ", sigma=", sigma) |> factor(levels = paste0("nu=", nu, ", mu=", round(mu_vals, 2), ", sigma=", sigma)) # フレーム切替用のラベル ) |> # 統計量を計算 dplyr::mutate( sd_minus = mean - sd, sd_plus = mean + sd ) |> # 期待値±標準偏差を計算 dplyr::select(!c(nu, mu, sd)) |> # 不要な列を削除 tidyr::pivot_longer( cols = !parameter, names_to = "type", values_to = "statistic" ) |> # 統計量の列をまとめる dplyr::mutate( type = stringr::str_replace(type, pattern = "sd_.*", replacement = "sd")) # 期待値±標準偏差のカテゴリを統一 anime_stat_df
## # A tibble: 183 × 3 ## parameter type statistic ## <fct> <chr> <dbl> ## 1 nu=5, mu=-3, sigma=0.5 mean -3 ## 2 nu=5, mu=-3, sigma=0.5 sd -3.65 ## 3 nu=5, mu=-3, sigma=0.5 sd -2.35 ## 4 nu=5, mu=-2.9, sigma=0.5 mean -2.9 ## 5 nu=5, mu=-2.9, sigma=0.5 sd -3.55 ## 6 nu=5, mu=-2.9, sigma=0.5 sd -2.25 ## 7 nu=5, mu=-2.8, sigma=0.5 mean -2.8 ## 8 nu=5, mu=-2.8, sigma=0.5 sd -3.45 ## 9 nu=5, mu=-2.8, sigma=0.5 sd -2.15 ## 10 nu=5, mu=-2.7, sigma=0.5 mean -2.7 ## # … with 173 more rows
「自由度の影響」のnu_valsをnu置き換え、mu_valsを追加するように処理します。
「自由度の影響」のコードで作図できます。フレーム数はmu_valsの要素数です。
# gif画像を作成 gganimate::animate(anime_dens_graph, nframes = length(mu_vals), fps = 10, width = 800, height = 600)
スケールパラメータの影響
$\sigma$の値を変化させ、$\nu, \mu$を固定して、anime_dens_dfを作成します。
・作図コード(クリックで展開)
パラメータごとに統計量を計算します。
# パラメータごとに統計量を計算 anime_stat_df <- tibble::tibble( nu = nu, sigma = sigma_vals, mean = dplyr::if_else( nu > 1, true = mu, false = as.numeric(NA) ), # 期待値 sd = dplyr::if_else( nu > 2, true = sqrt(sigma^2 * nu / (nu - 2)), false = as.numeric(NA) ), # 標準偏差 parameter = paste0("nu=", nu, ", mu=", mu, ", sigma=", round(sigma_vals, 2)) |> factor(levels = paste0("nu=", nu, ", mu=", mu, ", sigma=", round(sigma_vals, 2))) # フレーム切替用のラベル ) |> # 統計量を計算 dplyr::mutate( sd_minus = mean - sd, sd_plus = mean + sd ) |> # 期待値±標準偏差を計算 dplyr::select(!c(nu, sigma, sd)) |> # 不要な列を削除 tidyr::pivot_longer( cols = !parameter, names_to = "type", values_to = "statistic" ) |> # 統計量の列をまとめる dplyr::mutate( type = stringr::str_replace(type, pattern = "sd_.*", replacement = "sd")) # 期待値±標準偏差のカテゴリを統一 anime_stat_df
## # A tibble: 138 × 3 ## parameter type statistic ## <fct> <chr> <dbl> ## 1 nu=5, mu=2, sigma=0.5 mean 2 ## 2 nu=5, mu=2, sigma=0.5 sd 1.35 ## 3 nu=5, mu=2, sigma=0.5 sd 2.65 ## 4 nu=5, mu=2, sigma=0.6 mean 2 ## 5 nu=5, mu=2, sigma=0.6 sd 1.23 ## 6 nu=5, mu=2, sigma=0.6 sd 2.77 ## 7 nu=5, mu=2, sigma=0.7 mean 2 ## 8 nu=5, mu=2, sigma=0.7 sd 1.10 ## 9 nu=5, mu=2, sigma=0.7 sd 2.90 ## 10 nu=5, mu=2, sigma=0.8 mean 2 ## # … with 128 more rows
「位置パラメータの影響」のmu_valsとsigma_valsを置き換えるように処理します。
「自由度の影響」のコードで作図できます。フレーム数はsigma_valsの要素数です。
# gif画像を作成 gganimate::animate(anime_dens_graph, nframes = length(sigma_vals), fps = 10, width = 800, height = 600)
逆スケールパラメータの影響
$\sigma$の代わりに$\lambda$の値を変化させ、$\nu, \mu$を固定して、anime_dens_dfを作成します。
・作図コード(クリックで展開)
パラメータごとに統計量を計算します。
# パラメータごとに統計量を計算 anime_stat_df <- tibble::tibble( nu = nu, lambda = lambda_vals, mean = dplyr::if_else( nu > 1, true = mu, false = as.numeric(NA) ), # 期待値 sd = dplyr::if_else( nu > 2, true = nu / (nu - 2) / lambda, false = as.numeric(NA) ), # 標準偏差 parameter = paste0("nu=", nu, ", mu=", mu, ", lambda=", round(lambda_vals, 2)) |> factor(levels = paste0("nu=", nu, ", mu=", mu, ", lambda=", round(lambda_vals, 2))) # フレーム切替用のラベル ) |> # 統計量を計算 dplyr::mutate( sd_minus = mean - sd, sd_plus = mean + sd ) |> # 期待値±標準偏差を計算 dplyr::select(!c(nu, lambda, sd)) |> # 不要な列を削除 tidyr::pivot_longer( cols = !parameter, names_to = "type", values_to = "statistic" ) |> # 統計量の列をまとめる dplyr::mutate( type = stringr::str_replace(type, pattern = "sd_.*", replacement = "sd")) # 期待値±標準偏差のカテゴリを統一 anime_stat_df
## # A tibble: 300 × 3 ## parameter type statistic ## <fct> <chr> <dbl> ## 1 nu=5, mu=2, lambda=0.1 mean 2 ## 2 nu=5, mu=2, lambda=0.1 sd -14.7 ## 3 nu=5, mu=2, lambda=0.1 sd 18.7 ## 4 nu=5, mu=2, lambda=0.2 mean 2 ## 5 nu=5, mu=2, lambda=0.2 sd -6.33 ## 6 nu=5, mu=2, lambda=0.2 sd 10.3 ## 7 nu=5, mu=2, lambda=0.3 mean 2 ## 8 nu=5, mu=2, lambda=0.3 sd -3.56 ## 9 nu=5, mu=2, lambda=0.3 sd 7.56 ## 10 nu=5, mu=2, lambda=0.4 mean 2 ## # … with 290 more rows
「スケールパラメータの影響」のsigma_valsとlambda_valsを置き換えるように処理します。分散の計算も$\lambda$を使った式に変更します。
「自由度の影響」のコードで作図できます。フレーム数はlambda_valsの要素数です。
# gif画像を作成 gganimate::animate(anime_dens_graph, nframes = length(lambda_vals), fps = 10, width = 800, height = 600)
この記事では、1次元スチューデントのt分布の作図を確認しました。
参考文献
- C.M.ビショップ著,元田 浩・他訳『パターン認識と機械学習 上』,丸善出版,2012年.
おわりに
PRMLを参考にここまでの数記事を書いたので、ひょっとして久しく進んでいないPRMLの進捗と言えるのかな。
【次の内容】
つづく
