はじめに
機械学習で登場する確率分布について色々な角度から理解したいシリーズです。
この記事では、R言語でベータ分布のグラフを作成します。
【前の内容】
【他の記事一覧】
【この記事の内容】
ベータ分布の作図
ベータ分布(Beta Distribution)のグラフを作成します。ベータ分布については「1.2.3：ベータ分布【『トピックモデル』の勉強ノート】 - からっぽのしょこ」を参照してください。
利用するパッケージを読み込みます。
# 利用するパッケージ library(tidyverse) library(gganimate)
この記事では、基本的にパッケージ名::関数名()の記法を使うので、パッケージを読み込む必要はありません。ただし、作図コードがごちゃごちゃしないようにパッケージ名を省略しているため、ggplot2は読み込む必要があります。また、基本的にベースパイプ(ネイティブパイプ)演算子|>を使いますが、パイプ演算子%>%ででないと処理できない部分があるため、magrittrも読み込む必要があります。
分布の変化をアニメーション(gif画像)で確認するのにgganimateパッケージを利用します。不要であれば省略してください。
定義式の確認
まずは、ベータ分布の定義式を確認します。
ベータ分布は、次の式で定義されます。
ここで、$\alpha, \beta$はパラメータで、$\alpha > 0, \beta > 0$の値を満たす必要があります。確率変数の実現値$\phi$は、$0 < \phi < 1$となります。
ベータ分布は、0から1の値を生成することから、ベルヌーイ分布と二項分布のパラメータ(成功確率)の生成分布や事前分布として利用されます。
ベータ分布の期待値・分散・最頻値は、次の式で計算できます。詳しくは「1.2.3：ベータ分布【『トピックモデル』の勉強ノート】 - からっぽのしょこ」を参照してください。
ベータ分布の歪度と尖度は、次の式で計算できます。
これらの計算を行いグラフを作成します。
グラフの作成
ggplot2パッケージを利用して、ベータ分布のグラフを作成します。ベータ分布の確率や統計量の計算については「【R】ベータ分布の計算 - からっぽのしょこ」を参照してください。
ベータ分布のパラメータ$\alpha, \beta$を設定します。
# パラメータを指定 alpha <- 5 beta <- 2
$\alpha > 0, \beta > 0$の値を指定します。
ベータ分布の確率変数$\phi$がとり得る値を作成します。
# phiがとり得る値を作成 phi_vals <- seq(from = 0, to = 1, length.out = 501)
$0 \leq \phi \leq 1$の値を作成してphi_valsとします。$\phi$がとり得る範囲外の値の場合は、確率密度が0になります。
$\phi$の値ごとに確率密度を計算します。
# ベータ分布を計算 dens_df <- tidyr::tibble( phi = phi_vals, # 確率変数 density = dbeta(x = phi_vals, shape1 = alpha, shape2 = beta) # 確率密度 ) dens_df
## # A tibble: 501 × 2 ## phi density ## <dbl> <dbl> ## 1 0 0 ## 2 0.002 4.79e-10 ## 3 0.004 7.65e- 9 ## 4 0.006 3.86e- 8 ## 5 0.008 1.22e- 7 ## 6 0.01 2.97e- 7 ## 7 0.012 6.15e- 7 ## 8 0.014 1.14e- 6 ## 9 0.016 1.93e- 6 ## 10 0.018 3.09e- 6 ## # … with 491 more rows
ベータ分布の確率密度は、dbeta()で計算できます。確率変数の引数xにphi_vals、パラメータの引数shape1, shape2にalpha, betaを指定します。
phi_valsと、phi_valsの各要素に対応する確率密度をデータフレームに格納します。
ベータ分布のグラフを作成します。
# ベータ分布を作図 ggplot(data = dens_df, mapping = aes(x = phi, y = density)) + # データ geom_line(color = "#00A968", size = 1) + # 折れ線グラフ labs(title = "Beta Distribution", subtitle = paste0("alpha=", alpha, ", beta=", beta), # (文字列表記用) #subtitle = parse(text = paste0("list(alpha==", alpha, ", beta==", beta, ")")), # (数式表記用) x = expression(phi), y = "density") # ラベル
ギリシャ文字などの記号を使った数式を表示する場合は、expression()の記法を使います。等号は"=="、複数の(数式上の)変数を並べるには"list(変数1, 変数2)"とします。(プログラム上の)変数の値を使う場合は、parse()のtext引数に指定します。
この分布に統計量の情報を重ねて表示します。
# 補助線用の統計量を計算 E_phi <- alpha / (alpha + beta) s_phi <- sqrt(alpha * beta / (alpha + beta)^2 / (alpha + beta + 1)) mode_phi <- (alpha - 1) / (alpha + beta - 2) # 統計量を重ねたベータ分布のグラフを作成:線のみ ggplot(data = dens_df, mapping = aes(x = phi, y = density)) + # データ geom_line(color = "#00A968", size = 1) + # 折れ線グラフ geom_vline(xintercept = E_phi, color = "blue", size = 1, linetype = "dashed") + # 期待値 geom_vline(xintercept = E_phi-s_phi, color = "orange", size = 1, linetype = "dotted") + # 期待値 - 標準偏差 geom_vline(xintercept = E_phi+s_phi, color = "orange", size = 1, linetype = "dotted") + # 期待値 + 標準偏差 geom_vline(xintercept = mode_phi, color = "chocolate", size = 1, linetype = "dashed") + # 最頻値 labs(title = "Beta Distribution", subtitle = parse(text = paste0("list(alpha==", alpha, ", beta==", beta, ")")), x = expression(phi), y = "density") # ラベル
期待値・標準偏差(分散の平方根)・最頻値を計算して、それぞれgeom_vline()で垂直線を引きます。標準偏差については期待値から足し引きした値を使います。
凡例を表示する場合は、垂直線を引く値をデータフレームにまとめます。
# 統計量を格納 stat_df <- tibble::tibble( statistic = c(E_phi, E_phi-s_phi, E_phi+s_phi, mode_phi), # 統計量 type = c("mean", "sd", "sd", "mode") # 色分け用ラベル ) stat_df
## # A tibble: 4 × 2 ## statistic type ## <dbl> <chr> ## 1 0.714 mean ## 2 0.555 sd ## 3 0.874 sd ## 4 0.8 mode
各統計量を区別するための文字列も格納します。期待値と標準偏差の差・和は同じタイプとします。
線の色や種類、凡例テキストを指定する場合は、設定用の名前付きベクトルを作成します。
# 凡例用の設定を作成:(数式表示用) color_vec <- c(mean = "blue", sd = "orange", mode = "chocolate") linetype_vec <- c(mean = "dashed", sd = "dotted", mode = "dashed") label_vec <- c(mean = expression(E(x)), sd = expression(E(x) %+-% sqrt(V(x))), mode = expression(mode(x))) color_vec; linetype_vec; label_vec
## mean sd mode ## "blue" "orange" "chocolate" ## mean sd mode ## "dashed" "dotted" "dashed" ## expression(mean = E(x), sd = E(x) %+-% sqrt(V(x)), mode = mode(x))
各要素(設定)の名前をtype列の文字列と対応させて、設定(線の色・種類と数式用のテキスト)を指定します。
凡例付きのグラフを作成します。
# 統計量を重ねたベータ分布のグラフを作成:凡例付き ggplot() + geom_line(data = dens_df, mapping = aes(x = phi, 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") + # 線の種類:(線指定と数式表示用) theme(legend.text.align = 0) + # 図の体裁:凡例 labs(title = "Beta Distribution", subtitle = parse(text = paste0("list(alpha==", alpha, ", beta==", beta, ")")), x = expression(phi), y = "density") # ラベル
scale_color_manual()のvalues引数に線の色、scale_linetype_manual()のvalues引数に線の種類を指定します。凡例テキストを指定する場合は、それぞれのlabels引数に指定します。names引数は凡例ラベルです。
凡例テキストは、theme()のlegend.text.align引数に0を指定すると左寄せ、デフォルト(1)だと右寄せになります。
$\alpha \neq \beta$のとき、ベータ分布は非対称な形状になるので、期待値(オレンジ色の破線)と最頻値(茶色の破線)
が一致しません。
ベータ分布のグラフを描画できました。以降は、ここまでの作図処理を用いて、パラメータの影響を確認していきます。
パラメータと分布の関係を並べて比較
複数のパラメータのグラフを比較することで、パラメータの値と分布の形状の関係を確認します。
αの影響
複数の$\alpha$を指定し、$\beta$を固定して、それぞれ分布を計算します。
・作図コード(クリックで展開)
# パラメータとして利用する値を指定 alpha_vals <- c(0.1, 0.5, 1, 2, 4, 8) # 固定するパラメータを指定 beta <- 2 # phiがとり得る値を作成 phi_vals <- seq(from = 0, to = 1, length.out = 501) # パラメータごとにベータ分布を計算 res_dens_df <- tidyr::expand_grid( phi = phi_vals, alpha = alpha_vals ) |> # 全ての組み合わせを作成 dplyr::arrange(alpha, phi) |> # パラメータごとに並べ替え dplyr::mutate( dens = dbeta(x = phi, shape1 = alpha, shape2 = beta), parameter = paste0("alpha=", alpha, ", beta=", beta) |> factor(levels = paste0("alpha=", sort(alpha_vals), ", beta=", beta)) # 色分け用ラベル ) # 確率密度を計算 res_dens_df
## # A tibble: 3,006 × 4 ## phi alpha dens parameter ## <dbl> <dbl> <dbl> <fct> ## 1 0 0.1 Inf alpha=0.1, beta=2 ## 2 0.002 0.1 29.5 alpha=0.1, beta=2 ## 3 0.004 0.1 15.8 alpha=0.1, beta=2 ## 4 0.006 0.1 10.9 alpha=0.1, beta=2 ## 5 0.008 0.1 8.42 alpha=0.1, beta=2 ## 6 0.01 0.1 6.87 alpha=0.1, beta=2 ## 7 0.012 0.1 5.82 alpha=0.1, beta=2 ## 8 0.014 0.1 5.06 alpha=0.1, beta=2 ## 9 0.016 0.1 4.47 alpha=0.1, beta=2 ## 10 0.018 0.1 4.02 alpha=0.1, beta=2 ## # … with 2,996 more rows
確率変数phi_valsとパラメータalpha_valsそれぞれの要素の全ての組み合わせをexpand_grid()で作成します。
組み合わせごとに確率密度を計算して、パラメータごとにラベルを作成します。ラベルは、グラフの設定や凡例テキストとして使います。文字列型だと文字列の基準で順序が決まるので、因子型にしてパラメータに応じたレベル(順序)を設定します。
凡例を数式で表示する場合は、expression()に対応した記法に変換します。
# 凡例用のラベルを作成:(数式表示用) label_vec <- res_dens_df[["parameter"]] |> unique() |> # 重複を除去 stringr::str_replace_all(pattern = "=", replacement = "==") %>% # 等号表示用の記法に変換 paste0("list(", ., ")") |> # カンマ表示用の記法に変換 parse(text = _) # expression関数化 names(label_vec) <- unique(res_dens_df[["parameter"]]) # ggplotに指定する文字列に対応する名前付きベクトルに変換 label_vec[1]
## expression(`alpha=0.1, beta=2` = list(alpha == 0.1, beta == 2))
変換後の文字列のベクトルに対してnames()を使って、各要素の名前として元の文字列を設定します。
res_dens_dfのparameter列の要素alpha=0.1, beta=2と、変換後の文字列list(alpha == 0.1, beta == 2)が対応しています。
パラメータごとにベータ分布を作図します。
# パラメータごとにベータ分布を作図 ggplot(data = res_dens_df, mapping = aes(x = phi, y = dens, color = parameter)) + # データ geom_line(size = 1, show.legend = TRUE) + # 折れ線グラフ scale_color_hue(labels = label_vec) + # 線の色:(数式表示用) theme(legend.text.align = 0) + # 図の体裁:凡例 coord_cartesian(ylim = c(0, 5)) + # 軸の表示範囲 labs(title = "Beta Distribution", x = expression(phi), y = "density") # ラベル
βの影響
続いて、$\alpha$を固定し、複数の$\beta$を指定して、それぞれ分布を計算します。
・作図コード(クリックで展開)
# 固定するパラメータを指定 alpha <- 2 # パラメータとして利用する値を指定 beta_vals <- c(0.1, 0.5, 1, 2, 4, 8) # パラメータごとにベータ分布を計算 res_dens_df <- tidyr::expand_grid( phi = phi_vals, beta = beta_vals ) |> # 全ての組み合わせを作成 dplyr::arrange(beta, phi) |> # パラメータごとに並べ替え dplyr::mutate( dens = dbeta(x = phi, shape1 = alpha, shape2 = beta), parameter = paste0("alpha=", alpha, ", beta=", beta) |> factor(levels = paste0("alpha=", alpha, ", beta=", sort(beta_vals))) # 色分け用ラベル ) # 確率密度を計算 res_dens_df
## # A tibble: 3,006 × 4 ## phi beta dens parameter ## <dbl> <dbl> <dbl> <fct> ## 1 0 0.1 0 alpha=2, beta=0.1 ## 2 0.002 0.1 0.000220 alpha=2, beta=0.1 ## 3 0.004 0.1 0.000442 alpha=2, beta=0.1 ## 4 0.006 0.1 0.000664 alpha=2, beta=0.1 ## 5 0.008 0.1 0.000886 alpha=2, beta=0.1 ## 6 0.01 0.1 0.00111 alpha=2, beta=0.1 ## 7 0.012 0.1 0.00133 alpha=2, beta=0.1 ## 8 0.014 0.1 0.00156 alpha=2, beta=0.1 ## 9 0.016 0.1 0.00179 alpha=2, beta=0.1 ## 10 0.018 0.1 0.00201 alpha=2, beta=0.1 ## # … with 2,996 more rows
phi_valsとbeta_valsの全ての組み合わせを作成して、同様に処理します。
作図については同じコードで処理できます。
αとβの影響
今度は、$\alpha, \beta$の組み合わせを指定して、それぞれ分布を計算します。
・作図コード(クリックで展開)
# パラメータとして利用する値を指定 alpha_vals <- c(1, 0.5, 2, 2, 0.9, 0.8) beta_vals <- c(1, 0.5, 2, 4, 0.7, 1.2) # パラメータごとにベータ分布を計算 res_dens_df <- tidyr::expand_grid( phi = phi_vals, i = 1:length(alpha_vals) # パラメータ番号 ) |> # 全ての組み合わせを作成 dplyr::arrange(i, phi) |> # パラメータごとに並べ替え dplyr::mutate( alpha = alpha_vals[i], beta = beta_vals[i] ) |> # パラメータ列を追加 dplyr::mutate( dens = dbeta(x = phi, shape1 = alpha, shape2 = beta), parameter = paste0("alpha=", alpha, ", beta=", beta) |> factor() # 色分け用ラベル ) # 確率密度を計算 res_dens_df
## # A tibble: 3,006 × 6 ## phi i alpha beta dens parameter ## <dbl> <int> <dbl> <dbl> <dbl> <fct> ## 1 0 1 1 1 1 alpha=1, beta=1 ## 2 0.002 1 1 1 1 alpha=1, beta=1 ## 3 0.004 1 1 1 1 alpha=1, beta=1 ## 4 0.006 1 1 1 1 alpha=1, beta=1 ## 5 0.008 1 1 1 1 alpha=1, beta=1 ## 6 0.01 1 1 1 1 alpha=1, beta=1 ## 7 0.012 1 1 1 1 alpha=1, beta=1 ## 8 0.014 1 1 1 1 alpha=1, beta=1 ## 9 0.016 1 1 1 1 alpha=1, beta=1 ## 10 0.018 1 1 1 1 alpha=1, beta=1 ## # … with 2,996 more rows
同じ要素数となるようにalpha_vals, beta_valsに値を指定します。同じインデックスの値を使って分布を計算します。
phi_valsとパラメータのインデックスに対応する値の全ての組み合わせを作成して、alpha_vals, beta_valsから値を取り出し、それぞれ確率密度を計算します。
作図については同じコードで処理できます。
$\alpha = \beta$のとき、期待値と最頻値が$\phi = 0.5$になりますが、値によって分布の形状が異なります。また、$\alpha < 1, \beta < 1$のとき$\phi = 0, 1$で発散していることから、$\alpha > 1, \beta > 1$でのみ最頻値が求められるのを確認できます。
$\alpha$と$\beta$を比較して、$\alpha$が大きいほど$\phi$が1付近の確率密度が大きくなり、$\beta$が大きいほど$\phi$が0付近の確率密度が大きくなります。ただし、$\alpha > 1$のとき$\phi = 0$の確率密度が0、$\beta > 1$のとき$\phi = 1$の確率密度が0になります。
パラメータと分布の関係をアニメーションで可視化
前節では、複数のパラメータのグラフを並べて比較しました。次は、パラメータの値を少しずつ変化させて、分布の形状の変化をアニメーションで確認します。
αの影響
$\alpha$の値を変化させ、$\beta$を固定して、それぞれ分布を計算します。
・作図コード(クリックで展開)
# パラメータとして利用する値を指定 alpha_vals <- seq(from = 0.1, to = 15, by = 0.1) #length(alpha_vals) # フレーム数 # 固定するパラメータを指定 beta <- 2 # phiがとり得る値を作成 phi_vals <- seq(from = 0, to = 1, length.out = 501) # パラメータごとにベータ分布を計算 anime_dens_df <- tidyr::expand_grid( phi = phi_vals, alpha = alpha_vals ) |> # 全ての組み合わせを作成 dplyr::arrange(alpha, phi) |> # パラメータごとに並べ替え dplyr::mutate( dens = dbeta(x = phi, shape1 = alpha, shape2 = beta), parameter = paste0("alpha=", alpha, ", beta=", beta) |> factor(levels = paste0("alpha=", sort(alpha_vals), ", beta=", beta)) # フレーム切替用ラベル ) # 確率密度を計算 anime_dens_df
## # A tibble: 75,150 × 4 ## phi alpha dens parameter ## <dbl> <dbl> <dbl> <fct> ## 1 0 0.1 Inf alpha=0.1, beta=2 ## 2 0.002 0.1 29.5 alpha=0.1, beta=2 ## 3 0.004 0.1 15.8 alpha=0.1, beta=2 ## 4 0.006 0.1 10.9 alpha=0.1, beta=2 ## 5 0.008 0.1 8.42 alpha=0.1, beta=2 ## 6 0.01 0.1 6.87 alpha=0.1, beta=2 ## 7 0.012 0.1 5.82 alpha=0.1, beta=2 ## 8 0.014 0.1 5.06 alpha=0.1, beta=2 ## 9 0.016 0.1 4.47 alpha=0.1, beta=2 ## 10 0.018 0.1 4.02 alpha=0.1, beta=2 ## # … with 75,140 more rows
値の間隔が一定になるようにalpha_valsを作成します。パラメータごとにフレームを切り替えるので、alpha_valsの要素数がアニメーションのフレーム数になります。
データフレームの変数名以外は「並べて比較」のときと同じ処理です。アニメーションの作図では、parameter列をフレーム切替用のラベルとして使います。
ベータ分布のアニメーション(gif画像)を作成します。
# ベータ分布のアニメーションを作図 anime_dens_graph <- ggplot(data = anime_dens_df, mapping = aes(x = phi, y = dens)) + # データ geom_line(color = "#00A968", size = 1) + # 折れ線グラフ gganimate::transition_manual(parameter) + # フレーム coord_cartesian(ylim = c(0, 5)) + # 軸の表示範囲 labs(title = "Beta Distribution", subtitle = "{current_frame}", x = expression(phi), y = "density") # ラベル # gif画像を作成 gganimate::animate(anime_dens_graph, nframes = length(alpha_vals), fps = 10, width = 800, height = 600)
transition_manual()にフレームの順序を表す列を指定します。この例では、因子型のラベルのレベルの順に描画されます。
animate()のフレーム数の引数nframesにパラメータ数、フレームレートの引数fpsに1秒当たりのフレーム数を指定します。fps引数の値が大きいほどフレームが早く切り替わります。ただし、値が大きいと意図通りに動作しません。
βの影響
続いて、$\alpha$を固定し、$\beta$の値を変化させて、それぞれ分布を計算します。
・作図コード(クリックで展開)
# 固定するパラメータを指定 alpha <- 2 # パラメータとして利用する値を指定 beta_vals <- seq(from = 0.1, to = 15, by = 0.1) #length(beta_vals) # フレーム数 # パラメータごとにベータ分布を計算 anime_dens_df <- tidyr::expand_grid( phi = phi_vals, beta = beta_vals ) |> # 全ての組み合わせを作成 dplyr::arrange(beta, phi) |> # パラメータごとに並べ替え dplyr::mutate( dens = dbeta(x = phi, shape1 = alpha, shape2 = beta), parameter = paste0("alpha=", alpha, ", beta=", beta) |> factor(levels = paste0("alpha=", alpha, ", beta=", sort(beta_vals))) # フレーム切替用ラベル ) # 確率密度を計算 anime_dens_df
## # A tibble: 75,150 × 4 ## phi beta dens parameter ## <dbl> <dbl> <dbl> <fct> ## 1 0 0.1 0 alpha=2, beta=0.1 ## 2 0.002 0.1 0.000220 alpha=2, beta=0.1 ## 3 0.004 0.1 0.000442 alpha=2, beta=0.1 ## 4 0.006 0.1 0.000664 alpha=2, beta=0.1 ## 5 0.008 0.1 0.000886 alpha=2, beta=0.1 ## 6 0.01 0.1 0.00111 alpha=2, beta=0.1 ## 7 0.012 0.1 0.00133 alpha=2, beta=0.1 ## 8 0.014 0.1 0.00156 alpha=2, beta=0.1 ## 9 0.016 0.1 0.00179 alpha=2, beta=0.1 ## 10 0.018 0.1 0.00201 alpha=2, beta=0.1 ## # … with 75,140 more rows
「並べて比較」や「αの影響」と同様に処理します。
「αの影響」のコードで作図できます。フレーム数はbeta_valsの要素数です。
gganimate::animate(anime_dens_graph, nframes = length(beta_vals), fps = 10, width = 800, height = 600)
2パラメータの影響をアニメーションで可視化
前節では、1つのパラメータを固定して、もう1つのパラメータを変化させました。次は、固定するパラメータを複数個指定し、並べてアニメーションを作成します。
αの影響
複数の$\alpha$を指定し、$\beta$の値を変化させて、それぞれ分布を計算します。
・作図コード(クリックで展開)
# 比較する値を指定 alpha_vals <- c(0.1, 0.25, 0.5, 1, 1.5, 2, 5, 10, 15) # 変化する値を指定 beta_vals <- seq(from = 0.5, to = 10, by = 0.5) #length(beta_vals) # フレーム数 # phiがとり得る値を作成 phi_vals <- seq(from = 0, to = 1, length.out = 501) # パラメータごとにベータ分布を計算 anime_dens_df <- tidyr::expand_grid( phi = phi_vals, alpha = alpha_vals, beta = beta_vals ) |> # 全ての組み合わせを作成 dplyr::arrange(alpha, beta, phi) |> # パラメータごとに並べ替え dplyr::mutate( dens = dbeta(x = phi, shape1 = alpha, shape2 = beta) ) # 確率密度を計算 anime_dens_df
## # A tibble: 90,180 × 4 ## phi alpha beta dens ## <dbl> <dbl> <dbl> <dbl> ## 1 0 0.1 0.5 Inf ## 2 0.002 0.1 0.5 23.7 ## 3 0.004 0.1 0.5 12.7 ## 4 0.006 0.1 0.5 8.85 ## 5 0.008 0.1 0.5 6.84 ## 6 0.01 0.1 0.5 5.60 ## 7 0.012 0.1 0.5 4.76 ## 8 0.014 0.1 0.5 4.15 ## 9 0.016 0.1 0.5 3.68 ## 10 0.018 0.1 0.5 3.31 ## # … with 90,170 more rows
並べて比較する値をalpha_valsに指定します。
値の間隔が一定になるようにbeta_valsを作成します。
phi_vals・alpha_vals・beta_valsの全ての組み合わせを作成して、それぞれ確率密度を計算します。
パラメータごとに画面分割したアニメーション(gif画像)を作成します。
# ベータ分布のアニメーションを作図 anime_dens_graph <- ggplot(data = anime_dens_df, mapping = aes(x = phi, y = dens, color = as.factor(beta))) + # データ geom_line(show.legend = FALSE) + # 折れ線グラフ gganimate::transition_reveal(beta) + # フレーム facet_wrap(. ~ alpha, labeller = label_bquote(alpha==.(alpha))) + # グラフの分割 coord_cartesian(ylim = c(0, 6)) + # 軸の表示範囲 labs(title = "Beta Distribution", subtitle = "beta={frame_along}", x = expression(phi), y = "density") # ラベル # gif画像を作成 gganimate::animate(anime_dens_graph, nframes = length(beta_vals)+10, end_pause = 10, fps = 10, width = 1200, height = 900)
transition_reveal()を使ってフレームを切り替えると、過去フレームのグラフが表示され続けます。
facet_wrap()に列を指定すると、その列の値ごとにグラフを分割して描画できます。
animate()のend_pause引数を指定すると、最後のフレームでグラフが一時停止(最後のグラフを指定したフレーム数表示)します。
βの影響
続いて、$\alpha$の値を変化させ、複数の$\beta$を指定して、それぞれ分布を計算します。
・作図コード(クリックで展開)
# 変化する値を指定 alpha_vals <- seq(from = 0.5, to = 10, by = 0.5) #length(alpha_vals) # フレーム数 # 比較する値を指定 beta_vals <- c(0.1, 0.25, 0.5, 1, 1.5, 2, 5, 10, 15) # パラメータごとにベータ分布を計算 anime_dens_df <- tidyr::expand_grid( phi = phi_vals, alpha = alpha_vals, beta = beta_vals ) |> # 全ての組み合わせを作成 dplyr::arrange(beta, alpha, phi) |> # パラメータごとに並べ替え dplyr::mutate( dens = dbeta(x = phi, shape1 = alpha, shape2 = beta) ) # 確率密度を計算 anime_dens_df
## # A tibble: 90,180 × 4 ## phi alpha beta dens ## <dbl> <dbl> <dbl> <dbl> ## 1 0 0.5 0.1 Inf ## 2 0.002 0.5 0.1 1.98 ## 3 0.004 0.5 0.1 1.40 ## 4 0.006 0.5 0.1 1.15 ## 5 0.008 0.5 0.1 0.995 ## 6 0.01 0.5 0.1 0.891 ## 7 0.012 0.5 0.1 0.815 ## 8 0.014 0.5 0.1 0.756 ## 9 0.016 0.5 0.1 0.708 ## 10 0.018 0.5 0.1 0.669 ## # … with 90,170 more rows
値の間隔が一定になるようにalpha_valsを作成します。
並べて比較する値をbeta_valsに指定してします。
phi_vals・alpha_vals・beta_valsの全ての組み合わせを作成して、それそれ確率密度を計算します。
アニメーションを作成します。
# ベータ分布のアニメーションを作図 anime_dens_graph <- ggplot(data = anime_dens_df, mapping = aes(x = phi, y = dens, color = as.factor(alpha))) + # データ geom_line(show.legend = FALSE) + # 折れ線グラフ gganimate::transition_reveal(alpha) + # フレーム facet_wrap(. ~ beta, labeller = label_bquote(beta==.(beta))) + # グラフの分割 coord_cartesian(ylim = c(0, 6)) + # 軸の表示範囲 labs(title = "Beta Distribution", subtitle = "alpha={frame_along}", x = expression(phi), y = "density") # ラベル # gif画像を作成 gganimate::animate(anime_dens_graph, nframes = length(alpha_vals)+10, end_pause = 10, fps = 10, width = 1200, height = 900)
「αの影響」のalphaとbetaを置き換えて処理します。
パラメータと統計量の関係をアニメーションで可視化
ここまでは、パラメータと分布の関係を確認しました。次は、パラメータと統計量と歪度・尖度の関係をアニメーションで確認します。
αの影響
$\alpha$の値を変化させ、$\beta$を固定して、それぞれ歪度と尖度を計算し、フレーム切替用のラベルを作成します。
・作図コード(クリックで展開)
# パラメータとして利用する値を作成 alpha_vals <- seq(from = 0.1, to = 15, by = 0.1) # 固定するパラメータを指定 beta <- 5 # 歪度を計算 denom_vec <- 2 * (beta - alpha_vals) * sqrt(alpha_vals + beta + 1) numer_vec <- (alpha_vals + beta + 2) * sqrt(alpha_vals * beta) skewness_vec <- denom_vec / numer_vec # 尖度を計算 denom_vec <- 6 * ((alpha_vals - beta)^2 * (alpha_vals + beta + 1) - alpha_vals * beta * (alpha_vals + beta + 2)) numer_vec <- alpha_vals * beta * (alpha_vals + beta + 2) * (alpha_vals + beta + 3) kurtosis_vec <- denom_vec / numer_vec # ラベル用のテキストを作成 label_vec <- paste0( "alpha=", alpha_vals, ", beta=", beta, ", skewness=", round(skewness_vec, 3), ", kurtosis=", round(kurtosis_vec, 3) ) head(label_vec)
## [1] "alpha=0.1, beta=5, skewness=4.821, kurtosis=29.82" ## [2] "alpha=0.2, beta=5, skewness=3.32, kurtosis=13.785" ## [3] "alpha=0.3, beta=5, skewness=2.639, kurtosis=8.465" ## [4] "alpha=0.4, beta=5, skewness=2.224, kurtosis=5.822" ## [5] "alpha=0.5, beta=5, skewness=1.935, kurtosis=4.249" ## [6] "alpha=0.6, beta=5, skewness=1.717, kurtosis=3.212"
パラメータごとに歪度と尖度を計算して、フレーム切替用のラベルとして利用します。
パラメータごとに分布を計算します。
# phiがとり得る値を作成 phi_vals <- seq(from = 0, to = 1, length.out = 501) # パラメータごとにベータ分布を計算 anime_dens_df <- tidyr::expand_grid( phi = phi_vals, alpha = alpha_vals ) |> # 全ての組み合わせを作成 dplyr::arrange(alpha, phi) |> # パラメータごとに並べ替え dplyr::mutate( dens = dbeta(x = phi, shape1 = alpha, shape2 = beta), parameter = rep(label_vec, each = length(phi_vals)) |> factor(levels = label_vec) # フレーム切替用ラベル ) # 確率密度を計算 anime_dens_df
## # A tibble: 75,150 × 4 ## phi alpha dens parameter ## <dbl> <dbl> <dbl> <fct> ## 1 0 0.1 Inf alpha=0.1, beta=5, skewness=4.821, kurtosis=29.82 ## 2 0.002 0.1 32.6 alpha=0.1, beta=5, skewness=4.821, kurtosis=29.82 ## 3 0.004 0.1 17.3 alpha=0.1, beta=5, skewness=4.821, kurtosis=29.82 ## 4 0.006 0.1 11.9 alpha=0.1, beta=5, skewness=4.821, kurtosis=29.82 ## 5 0.008 0.1 9.14 alpha=0.1, beta=5, skewness=4.821, kurtosis=29.82 ## 6 0.01 0.1 7.41 alpha=0.1, beta=5, skewness=4.821, kurtosis=29.82 ## 7 0.012 0.1 6.24 alpha=0.1, beta=5, skewness=4.821, kurtosis=29.82 ## 8 0.014 0.1 5.39 alpha=0.1, beta=5, skewness=4.821, kurtosis=29.82 ## 9 0.016 0.1 4.74 alpha=0.1, beta=5, skewness=4.821, kurtosis=29.82 ## 10 0.018 0.1 4.23 alpha=0.1, beta=5, skewness=4.821, kurtosis=29.82 ## # … with 75,140 more rows
これまでと同様に処理します。フレーム切替用のラベル(parameter列)の値にはlabel_vecを使います。
パラメータごとに統計量を計算します。
# 統計量を格納 anime_stat_df <- tibble::tibble( mean = alpha_vals / (alpha_vals + beta), # 期待値 sd = sqrt(alpha_vals * beta / (alpha_vals + beta)^2 / (alpha_vals + beta + 1)), # 標準偏差 mode = (alpha_vals - 1) / (alpha_vals + beta - 2), # 最頻値 parameter = factor(label_vec, levels = label_vec) # フレーム切替用ラベル ) |> # 統計量を計算 dplyr::mutate( sd_minus = mean - sd, sd_plus = mean + sd ) |> # 期待値±標準偏差を計算 dplyr::select(!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: 600 × 3 ## parameter type statistic ## <fct> <chr> <dbl> ## 1 alpha=0.1, beta=5, skewness=4.821, kurtosis=29.82 mean 0.0196 ## 2 alpha=0.1, beta=5, skewness=4.821, kurtosis=29.82 mode -0.290 ## 3 alpha=0.1, beta=5, skewness=4.821, kurtosis=29.82 sd -0.0365 ## 4 alpha=0.1, beta=5, skewness=4.821, kurtosis=29.82 sd 0.0757 ## 5 alpha=0.2, beta=5, skewness=3.32, kurtosis=13.785 mean 0.0385 ## 6 alpha=0.2, beta=5, skewness=3.32, kurtosis=13.785 mode -0.25 ## 7 alpha=0.2, beta=5, skewness=3.32, kurtosis=13.785 sd -0.0388 ## 8 alpha=0.2, beta=5, skewness=3.32, kurtosis=13.785 sd 0.116 ## 9 alpha=0.3, beta=5, skewness=2.639, kurtosis=8.465 mean 0.0566 ## 10 alpha=0.3, beta=5, skewness=2.639, kurtosis=8.465 mode -0.212 ## # … with 590 more rows
期待値・標準偏差・最頻値を計算して、さらに期待値から標準偏差を引いた値と足した値を求めます。
期待値・標準偏差の和と差・最頻値の4つの列をpivot_longer()でまとめます。列名がラベルになるので、標準偏差の和と差のラベルを統一します。
線と凡例の設定用の名前付きベクトルを作成します。
# 凡例用の設定を作成:(数式表示用) color_vec <- c(mean = "blue", sd = "orange", mode = "chocolate") linetype_vec <- c(mean = "dashed", sd = "dotted", mode = "dashed") label_vec <- c(mean = expression(E(x)), sd = expression(E(x) %+-% sqrt(V(x))), mode = expression(mode(x)))
「グラフの作成」のときと同じ処理です。
統計量の情報を重ねたベータ分布のアニメーション(gif画像)を作成します。
# 統計量を重ねた分布のアニメーションを作図 anime_prob_graph <- ggplot() + # データ geom_line(data = anime_dens_df, mapping = aes(x = phi, y = dens), 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_linetype_manual(values = linetype_vec, labels = label_vec, name = "statistic") + # 線の種類:(線指定と数式表示用) scale_color_manual(values = color_vec, labels = label_vec, name = "statistic") + # 線の色:(色指定と数式表示用) theme(legend.text.align = 0) + # 図の体裁:凡例 coord_cartesian(xlim = c(0, 1), ylim = c(0, 5)) + # 軸の表示範囲 labs(title = "Beta Distribution", subtitle = "{current_frame}", x = expression(phi), y = "density") # ラベル # gif画像を作成 gganimate::animate(anime_prob_graph, nframes = length(alpha_vals), fps = 10, width = 1200, height = 600)
βの影響
続いて、$\alpha$を固定し、$\beta$の値を変化させて、それぞれ歪度と尖度を計算し、フレーム切替用のラベルを作成します。
・作図コード(クリックで展開)
# 固定するパラメータを指定 alpha <- 5 # パラメータとして利用する値を作成 beta_vals <- seq(from = 0.1, to = 15, by = 0.1) # 歪度を計算 denom_vec <- 2 * (beta_vals - alpha) * sqrt(alpha + beta_vals + 1) numer_vec <- (alpha+ beta_vals + 2) * sqrt(alpha * beta_vals) skewness_vec <- denom_vec / numer_vec # 尖度を計算 denom_vec <- 6 * ((alpha - beta_vals)^2 * (alpha + beta_vals + 1) - alpha * beta_vals * (alpha + beta_vals + 2)) numer_vec <- alpha * beta_vals * (alpha + beta_vals + 2) * (alpha + beta_vals + 3) kurtosis_vec <- denom_vec / numer_vec # ラベル用のテキストを作成 label_vec <- paste0( "alpha=", alpha, ", beta=", beta_vals, ", skewness=", round(skewness_vec, 3), ", kurtosis=", round(kurtosis_vec, 3) ) head(label_vec)
## [1] "alpha=5, beta=0.1, skewness=-4.821, kurtosis=29.82" ## [2] "alpha=5, beta=0.2, skewness=-3.32, kurtosis=13.785" ## [3] "alpha=5, beta=0.3, skewness=-2.639, kurtosis=8.465" ## [4] "alpha=5, beta=0.4, skewness=-2.224, kurtosis=5.822" ## [5] "alpha=5, beta=0.5, skewness=-1.935, kurtosis=4.249" ## [6] "alpha=5, beta=0.6, skewness=-1.717, kurtosis=3.212"
「αの影響」のalpha_valsをalpha、betaをbeta_valsに置き換えると処理できます。
パラメータごとに分布を計算します。
# パラメータごとにベータ分布を計算 anime_dens_df <- tidyr::expand_grid( phi = phi_vals, beta = beta_vals ) |> # 全ての組み合わせを作成 dplyr::arrange(beta, phi) |> # パラメータごとに並べ替え dplyr::mutate( dens = dbeta(x = phi, shape1 = alpha, shape2 = beta), parameter = rep(label_vec, each = length(phi_vals)) |> factor(levels = label_vec) # フレーム切替用ラベル ) # 確率密度を計算 anime_dens_df
## # A tibble: 75,150 × 4 ## phi beta dens parameter ## <dbl> <dbl> <dbl> <fct> ## 1 0 0.1 0 alpha=5, beta=0.1, skewness=-4.821, kurtosis=29.82 ## 2 0.002 0.1 1.96e-12 alpha=5, beta=0.1, skewness=-4.821, kurtosis=29.82 ## 3 0.004 0.1 3.14e-11 alpha=5, beta=0.1, skewness=-4.821, kurtosis=29.82 ## 4 0.006 0.1 1.59e-10 alpha=5, beta=0.1, skewness=-4.821, kurtosis=29.82 ## 5 0.008 0.1 5.05e-10 alpha=5, beta=0.1, skewness=-4.821, kurtosis=29.82 ## 6 0.01 0.1 1.23e- 9 alpha=5, beta=0.1, skewness=-4.821, kurtosis=29.82 ## 7 0.012 0.1 2.56e- 9 alpha=5, beta=0.1, skewness=-4.821, kurtosis=29.82 ## 8 0.014 0.1 4.76e- 9 alpha=5, beta=0.1, skewness=-4.821, kurtosis=29.82 ## 9 0.016 0.1 8.13e- 9 alpha=5, beta=0.1, skewness=-4.821, kurtosis=29.82 ## 10 0.018 0.1 1.31e- 8 alpha=5, beta=0.1, skewness=-4.821, kurtosis=29.82 ## # … with 75,140 more rows
label_vecを使って、これまでと同様に処理します。
パラメータごとに統計量を計算します。
# 統計量を格納 anime_stat_df <- tibble::tibble( mean = alpha / (alpha + beta_vals), # 期待値 sd = sqrt(alpha * beta_vals / (alpha + beta_vals)^2 / (alpha + beta_vals + 1)), # 標準偏差 mode = (alpha - 1) / (alpha + beta_vals - 2), # 最頻値 parameter = factor(label_vec, levels = label_vec) # フレーム切替用ラベル ) |> # 統計量を計算 dplyr::mutate( sd_minus = mean - sd, sd_plus = mean + sd ) |> # 期待値±標準偏差を計算 dplyr::select(!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: 600 × 3 ## parameter type statistic ## <fct> <chr> <dbl> ## 1 alpha=5, beta=0.1, skewness=-4.821, kurtosis=29.82 mean 0.980 ## 2 alpha=5, beta=0.1, skewness=-4.821, kurtosis=29.82 mode 1.29 ## 3 alpha=5, beta=0.1, skewness=-4.821, kurtosis=29.82 sd 0.924 ## 4 alpha=5, beta=0.1, skewness=-4.821, kurtosis=29.82 sd 1.04 ## 5 alpha=5, beta=0.2, skewness=-3.32, kurtosis=13.785 mean 0.962 ## 6 alpha=5, beta=0.2, skewness=-3.32, kurtosis=13.785 mode 1.25 ## 7 alpha=5, beta=0.2, skewness=-3.32, kurtosis=13.785 sd 0.884 ## 8 alpha=5, beta=0.2, skewness=-3.32, kurtosis=13.785 sd 1.04 ## 9 alpha=5, beta=0.3, skewness=-2.639, kurtosis=8.465 mean 0.943 ## 10 alpha=5, beta=0.3, skewness=-2.639, kurtosis=8.465 mode 1.21 ## # … with 590 more rows
「αの影響」のalpha_valsをalpha、betaをbeta_valsに置き換えると処理できます。
「αの影響」のコードで作図できます。フレーム数はbeta_valsの要素数です。
# gif画像を作成 gganimate::animate(anime_prob_graph, nframes = length(beta_vals), fps = 10, width = 800, height = 600)
この記事では、ベータ分布の作図を確認しました。次は、乱数を生成します。
参考文献
- 岩田具治『トピックモデル』(機械学習プロフェッショナルシリーズ)講談社,2015年.
おわりに
加筆修正の際に記事を分割しました。
【次の内容】
