はじめに
機械学習で登場する確率分布について色々な角度から理解したいシリーズです。
この記事では、R言語でガンマ分布のグラフを作成します。
【前の内容】
【他の記事一覧】
【この記事の内容】
ガンマ分布の作図
ガンマ分布(Gamma Distribution)のグラフを作成します。ガンマ分布についてはいつか書きます。
利用するパッケージを読み込みます。
# 利用パッケージ library(tidyverse) library(gganimate)
この記事では、基本的にパッケージ名::関数名()の記法を使うので、パッケージを読み込む必要はありません。ただし、作図コードがごちゃごちゃしないようにパッケージ名を省略しているため、ggplot2は読み込む必要があります。
magrittrパッケージのパイプ演算子%>%ではなく、ベースパイプ(ネイティブパイプ)演算子|>を使っています。%>%に置き換えても処理できます。
分布の変化をアニメーション(gif画像)で確認するのにgganimateパッケージを利用します。不要であれば省略してください。
定義式の確認
まずは、ガンマ分布の定義式を確認します。
ガンマ分布は、次の式で定義されます。
ここで、$a$は形状に関するパラメータ、$b$は尺度に関するパラメータです。
確率変数の実現値$\lambda$は、非負の実数$\lambda > 0$となります。パラメータ$a, b$は、$a > 0, b > 0$を満たす必要があります。
ガンマ分布は、非負の実数値を生成することから、ポアソン分布のパラメータや1次元ガウス分布の精度パラメータの生成分布や事前分布として利用されます。
ガンマ分布の期待値・分散・最頻値は、次の式で計算できます。詳しくはいつか書きます。
ガンマ分布の歪度と尖度は、次の式で計算できます。詳しくはいつか書きます。
これらの計算を行いグラフを作成します。
グラフの作成
ggplot2パッケージを利用して、ガンマ分布のグラフを作成します。ガンマ分布の確率や統計量の計算については「【R】ガンマ分布の計算 - からっぽのしょこ」を参照してください。
ガンマ分布のパラメータ$a, b$を設定します。
# パラメータを指定 a <- 5 b <- 2
$a > 0, b > 0$の値を指定します。
ガンマ分布の確率変数がとり得る値$\lambda$を作成します。
# lambdaの値を作成 lambda_vals <- seq(from = 0, to = a/b * 4, length.out = 251)
$\lambda > 0$の値を作成してlambda_valsとします。この例では、0から期待値の4倍を範囲とします。$\lambda$がとり得る範囲外の値の場合は、確率密度が0になります。
$\lambda$の値ごとの確率密度を計算します。
# ガンマ分布を計算 dens_df <- tidyr::tibble( lambda = lambda_vals, # 確率変数 density = dgamma(x = lambda_vals, shape = a, rate = b) # 確率密度 ) dens_df
## # A tibble: 251 × 2 ## lambda density ## <dbl> <dbl> ## 1 0 0 ## 2 0.04 0.00000315 ## 3 0.08 0.0000465 ## 4 0.12 0.000217 ## 5 0.16 0.000635 ## 6 0.2 0.00143 ## 7 0.24 0.00274 ## 8 0.28 0.00468 ## 9 0.32 0.00737 ## 10 0.36 0.0109 ## # … with 241 more rows
ガンマ分布の確率密度は、dgamma()で計算できます。確率変数の引数xにlambda_vals、形状の引数shapeにa、尺度の引数rateにbを指定します。
lambda_valsと、lambda_valsの各要素に対応する確率密度をデータフレームに格納します。
ガンマ分布のグラフを作成します。
# ガンマ分布のグラフを作成 ggplot(data = dens_df, mapping = aes(x = lambda, y = density)) + # データ geom_line(color = "#00A968", size = 1) + # 折れ線グラフ labs(title = "Gamma Distribution", subtitle = paste0("a=", a, ", b=", b), x = expression(lambda)) # ラベル
この分布に統計量の情報を重ねて表示します。
# 補助線用の統計量を計算 E_lambda <- a / b s_lambda <- sqrt(a / b^2) mode_lambda <- (a - 1) / b # 統計量を重ねたガンマ分布のグラフを作成:線のみ ggplot(data = dens_df, mapping = aes(x = lambda, y = density)) + # データ geom_line(color = "#00A968", size = 1) + # 分布 geom_vline(xintercept = E_lambda, color = "blue", size = 1, linetype = "dashed") + # 期待値 geom_vline(xintercept = E_lambda - s_lambda, color = "orange", size = 1, linetype = "dotted") + # 期待値 - 標準偏差 geom_vline(xintercept = E_lambda + s_lambda, color = "orange", size = 1, linetype = "dotted") + # 期待値 + 標準偏差 geom_vline(xintercept = mode_lambda, color = "chocolate", size = 1, linetype = "dashed") + # 最頻値 labs(title = "Gamma Distribution", subtitle = paste0("a=", a, ", b=", b), x = expression(lambda)) # ラベル
期待値・標準偏差(分散の平方根)・最頻値を計算して、それぞれgeom_vline()で垂直線を引きます。標準偏差については期待値から足し引きした値を使います。
凡例を表示する場合は、垂直線を引く値をデータフレームにまとめます。
# 統計量を格納 stat_df <- tibble::tibble( statistic = c(E_lambda, E_lambda-s_lambda, E_lambda+s_lambda, mode_lambda), # 統計量 type = c("mean", "sd", "sd", "mode") # 色分け用ラベル ) stat_df
## # A tibble: 4 × 2 ## statistic type ## <dbl> <chr> ## 1 2.5 mean ## 2 1.38 sd ## 3 3.62 sd ## 4 2 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 = lambda, 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 = "Gamma Distribution", subtitle = paste0("a=", a, ", b=", b), x = expression(lambda), y = "density") # ラベル
scale_color_manual()のvalues引数に線の色、scale_linetype_manual()のvalues引数に線の種類を指定します。凡例テキストを指定する場合は、それぞれのlabels引数に指定します。names引数は凡例ラベルです。
凡例テキストは、theme()のlegend.text.align引数に0を指定すると左寄せ、デフォルト(1)だと右寄せになります。
ガンマ分布は非対称な形状になるので、期待値(オレンジ色の破線)と最頻値(茶色の破線)が一致しません。
ガンマ分布のグラフを描画できました。以降は、ここまでの作図処理を用いて、パラメータの影響を確認していきます。
パラメータと分布の関係を並べて比較
複数のパラメータのグラフを比較することで、パラメータの値と分布の形状の関係を確認します。
aの影響
複数の$a$を指定し、$b$を固定して、それぞれ分布を計算します。
・作図コード(クリックで展開)
# パラメータとして利用する値を指定 a_vals <- c(0.1, 1, 2.5, 5, 10.5) # 固定するパラメータを指定 b <- 5 # lambdaの値を作成 lambda_vals <- seq(from = 0, to = 5, length.out = 501) # パラメータごとにガンマ分布を計算 res_dens_df <- tidyr::expand_grid( lambda = lambda_vals, a = a_vals ) |> # 全ての組み合わせを作成 dplyr::arrange(a, lambda) |> # パラメータごとに並べ替え dplyr::mutate( density = dgamma(x = lambda, shape = a, rate = b), parameter = paste0("a=", a, ", b=", b) |> factor(levels = paste0("a=", sort(a_vals), ", b=", b)) # 色分け用ラベル ) # 確率密度を計算 res_dens_df
## # A tibble: 2,505 × 4 ## lambda a density parameter ## <dbl> <dbl> <dbl> <fct> ## 1 0 0.1 Inf a=0.1, b=5 ## 2 0.01 0.1 7.41 a=0.1, b=5 ## 3 0.02 0.1 3.78 a=0.1, b=5 ## 4 0.03 0.1 2.49 a=0.1, b=5 ## 5 0.04 0.1 1.83 a=0.1, b=5 ## 6 0.05 0.1 1.43 a=0.1, b=5 ## 7 0.06 0.1 1.15 a=0.1, b=5 ## 8 0.07 0.1 0.953 a=0.1, b=5 ## 9 0.08 0.1 0.804 a=0.1, b=5 ## 10 0.09 0.1 0.688 a=0.1, b=5 ## # … with 2,495 more rows
確率変数lambda_valsとパラメータa_valsそれぞれの要素の全ての組み合わせをexpand_grid()で作成します。
組み合わせごとに確率密度を計算して、パラメータごとにラベルを作成します。ラベルは、グラフの設定や凡例テキストとして使います。文字列型だと文字列の基準で順序が決まるので、因子型にしてパラメータに応じたレベル(順序)を設定します。
パラメータごとにガンマ分布を作図します。
# パラメータごとにガンマ分布のグラフを作成 ggplot(data = res_dens_df, mapping = aes(x = lambda, y = density, color = parameter)) + # データ geom_line(size = 1) + # 折れ線グラフ coord_cartesian(ylim = c(0, 5)) + # 軸の表示範囲 theme(legend.text.align = 0) + # 図の体裁:凡例 labs(title = "Gamma Distribution", x = expression(lambda), y = "density") # ラベル
$a$が小さいほど$\lambda = 0$付近の確率密度が大きくなります。
bの影響
続いて、$a$を固定し、複数の$b$を指定して、それぞれ分布を計算します。
・作図コード(クリックで展開)
# 固定するパラメータを指定 a <- 5 # パラメータとして利用する値を指定 b_vals <- c(0.5, 1, 2.5, 5, 10.5, 25) # パラメータごとにガンマ分布を計算 res_dens_df <- tidyr::expand_grid( lambda = lambda_vals, b = b_vals ) |> # 全ての組み合わせを作成 dplyr::arrange(b, lambda) |> # パラメータごとに並べ替え dplyr::mutate( density = dgamma(x = lambda, shape = a, rate = b), parameter = paste0("a=", a, ", b=", b) |> factor(levels = paste0("a=", a, ", b=", sort(b_vals))) # 色分け用ラベル ) # 確率密度を計算 res_dens_df
## # A tibble: 3,006 × 4 ## lambda b density parameter ## <dbl> <dbl> <dbl> <fct> ## 1 0 0.5 0 a=5, b=0.5 ## 2 0.01 0.5 1.30e-11 a=5, b=0.5 ## 3 0.02 0.5 2.06e-10 a=5, b=0.5 ## 4 0.03 0.5 1.04e- 9 a=5, b=0.5 ## 5 0.04 0.5 3.27e- 9 a=5, b=0.5 ## 6 0.05 0.5 7.94e- 9 a=5, b=0.5 ## 7 0.06 0.5 1.64e- 8 a=5, b=0.5 ## 8 0.07 0.5 3.02e- 8 a=5, b=0.5 ## 9 0.08 0.5 5.12e- 8 a=5, b=0.5 ## 10 0.09 0.5 8.17e- 8 a=5, b=0.5 ## # … with 2,996 more rows
lambda_valsとb_valsの全ての組み合わせを作成して、同様に処理します。
作図については同じコードで処理できます。
$b$が大きいほど$\lambda = 0$付近の確率密度が大きくなります。
aとbの影響
今度は、$a, b$の組み合わせを指定して、それぞれ分布を計算します。
・作図コード(クリックで展開)
# パラメータとして利用する値を指定 a_vals <- c(0.5, 6, 1, 2.5, 5.5) b_vals <- c(0.1, 0.9, 1, 10, 5.5) # パラメータごとにガンマ分布を計算 res_dens_df <- tidyr::expand_grid( lambda = lambda_vals, i = 1:length(a_vals) # パラメータ番号 ) |> # 全ての組み合わせを作成 dplyr::arrange(i, lambda) |> # パラメータごとに並べ替え dplyr::mutate( a = a_vals[i], b = b_vals[i] ) |> # パラメータ列を追加 dplyr::mutate( density = dgamma(x = lambda, shape = a, rate = b), parameter = paste0("a=", a, ", b=", b) |> factor() # 色分け用ラベル ) # 確率密度を計算 res_dens_df
## # A tibble: 2,505 × 6 ## lambda i a b density parameter ## <dbl> <int> <dbl> <dbl> <dbl> <fct> ## 1 0 1 0.5 0.1 Inf a=0.5, b=0.1 ## 2 0.01 1 0.5 0.1 1.78 a=0.5, b=0.1 ## 3 0.02 1 0.5 0.1 1.26 a=0.5, b=0.1 ## 4 0.03 1 0.5 0.1 1.03 a=0.5, b=0.1 ## 5 0.04 1 0.5 0.1 0.889 a=0.5, b=0.1 ## 6 0.05 1 0.5 0.1 0.794 a=0.5, b=0.1 ## 7 0.06 1 0.5 0.1 0.724 a=0.5, b=0.1 ## 8 0.07 1 0.5 0.1 0.670 a=0.5, b=0.1 ## 9 0.08 1 0.5 0.1 0.626 a=0.5, b=0.1 ## 10 0.09 1 0.5 0.1 0.589 a=0.5, b=0.1 ## # … with 2,495 more rows
同じ要素数となるようにa_vals, b_valsに値を指定します。同じインデックスの値を使って分布を計算します。
lambda_valsとパラメータのインデックスに対応する値の全ての組み合わせを作成して、a_vals, b_valsから値を取り出し、それぞれ確率密度を計算します。
作図については同じコードで処理できます。
パラメータと分布の形状の関係をアニメーションで可視化
前節では、複数のパラメータのグラフを並べて比較しました。次は、パラメータの値を少しずつ変化させて、分布の形状の変化をアニメーションで確認します。
aの影響
$a$の値を変化させ、$b$を固定して、それぞれ分布を計算します。
・作図コード(クリックで展開)
# パラメータとして利用する値を指定 a_vals <- seq(from = 0.1, to = 10, by = 0.1) length(a_vals) # フレーム数 # 固定するパラメータを指定 b <- 2 # lambdaの値を作成 lambda_vals <- seq(from = 0, to = 5, length.out = 250) # パラメータごとにガンマ分布を計算 anime_dens_df <- tidyr::expand_grid( lambda = lambda_vals, a = a_vals ) |> # 全ての組み合わせを作成 dplyr::arrange(a, lambda) |> # パラメータごとに並べ替え dplyr::mutate( density = dgamma(x = lambda, shape = a, rate = b), parameter = paste0("a=", a, ", b=", b) |> factor(levels = paste0("a=", a_vals, ", b=", b)) # フレーム切替用ラベル ) # 確率密度を計算 anime_dens_df
## # A tibble: 25,000 × 4 ## lambda a density parameter ## <dbl> <dbl> <dbl> <fct> ## 1 0 0.1 Inf a=0.1, b=2 ## 2 0.0201 0.1 3.65 a=0.1, b=2 ## 3 0.0402 0.1 1.88 a=0.1, b=2 ## 4 0.0602 0.1 1.25 a=0.1, b=2 ## 5 0.0803 0.1 0.928 a=0.1, b=2 ## 6 0.100 0.1 0.729 a=0.1, b=2 ## 7 0.120 0.1 0.595 a=0.1, b=2 ## 8 0.141 0.1 0.497 a=0.1, b=2 ## 9 0.161 0.1 0.424 a=0.1, b=2 ## 10 0.181 0.1 0.366 a=0.1, b=2 ## # … with 24,990 more rows
値の間隔が一定になるようにa_valsを作成します。パラメータごとにフレームを切り替えるので、a_valsの要素数がアニメーションのフレーム数になります。
データフレームの変数名以外は「並べて比較」のときと同じ処理です。アニメーションの作図では、parameter列をフレーム切替用のラベルとして使います。
ガンマ分布のアニメーション(gif画像)を作成します。
# ガンマ分布のアニメーションを作図 anime_dens_graph <- ggplot(data = anime_dens_df, mapping = aes(x = lambda, y = density)) + # データ geom_line(color = "#00A968", size = 1) + # 折れ線グラフ gganimate::transition_manual(parameter) + # フレーム coord_cartesian(ylim = c(0, 5)) + # 軸の表示範囲 labs(title = "Gamma Distribution", subtitle = "{current_frame}", x = expression(lambda), y = "density") # ラベル # gif画像を作成 gganimate::animate(anime_dens_graph, nframes = length(a_vals), fps = 10, width = 800, height = 600)
transition_manual()にフレームの順序を表す列を指定します。この例では、因子型のラベルのレベルの順に描画されます。
animate()のフレーム数の引数nframesにパラメータ数、フレームレートの引数fpsに1秒当たりのフレーム数を指定します。fps引数の値が大きいほどフレームが早く切り替わります。ただし、値が大きいと意図通りに動作しません。
$a$が大きくなるに従って、$\lambda$が大きいほど確率密度が大きくなり(山が右に移動し)ます。これは、期待値と最頻値の計算式からも分かります。
bの影響
続いて、$a$を固定し、$b$の値を変化させて、それぞれ分布を計算します。
・作図コード(クリックで展開)
# 固定するパラメータを指定 a <- 2 # パラメータとして利用する値を指定 b_vals <- seq(from = 0.1, to = 10, by = 0.1) length(a_vals) # フレーム数 # 作図用のlambdaの点を作成 lambda_vals <- seq(from = 0, to = 5, length.out = 251) # パラメータごとにガンマ分布を計算 anime_dens_df <- tidyr::expand_grid( lambda = lambda_vals, b = b_vals ) |> # 全ての組み合わせを作成 dplyr::arrange(b, lambda) |> # パラメータごとに並べ替え dplyr::mutate( density = dgamma(x = lambda, shape = a, rate = b), parameter = paste0("a=", a, ", b=", b) |> factor(levels = paste0("a=", a, ", b=", b_vals)) # フレーム切替用ラベル ) # 確率密度を計算 anime_dens_df
## # A tibble: 25,100 × 4 ## lambda b density parameter ## <dbl> <dbl> <dbl> <fct> ## 1 0 0.1 0 a=2, b=0.1 ## 2 0.02 0.1 0.000200 a=2, b=0.1 ## 3 0.04 0.1 0.000398 a=2, b=0.1 ## 4 0.06 0.1 0.000596 a=2, b=0.1 ## 5 0.08 0.1 0.000794 a=2, b=0.1 ## 6 0.1 0.1 0.000990 a=2, b=0.1 ## 7 0.12 0.1 0.00119 a=2, b=0.1 ## 8 0.14 0.1 0.00138 a=2, b=0.1 ## 9 0.16 0.1 0.00157 a=2, b=0.1 ## 10 0.18 0.1 0.00177 a=2, b=0.1 ## # … with 25,090 more rows
「並べて比較」や「aの影響」と同様に処理します。
「aの影響」のコードで作図できます。フレーム数はb_valsの要素数です。
# gif画像を作成 gganimate::animate(anime_dens_graph, nframes = length(b_vals), fps = 10, width = 800, height = 600)
$b$が大きくなるに従って、分布が広がり(山が低くなり)ます。これは、分散の計算式からも分かります。
2パラメータの影響をアニメーションで可視化
前節では、1つのパラメータを固定して、もう1つのパラメータを変化させました。次は、固定するパラメータを複数個指定し、並べてアニメーションを作成します。
αの影響
複数の$a$を指定し、$b$の値を変化させて、それぞれ分布を計算します。
・作図コード(クリックで展開)
# 比較する値を指定 a_vals <- c(0.5, 1, 2, 4, 8, 16) # 変化する値を指定 b_vals <- seq(from = 0.5, to = 10, by = 0.5) length(b_vals) # フレーム数 # lambdaの値を作成 lambda_vals <- seq(from = 0, to = 5, length.out = 251) # パラメータごとにガンマ分布を計算 anime_dens_df <- tidyr::expand_grid( lambda = lambda_vals, a = a_vals, b = b_vals ) |> # 全ての組み合わせを作成 dplyr::arrange(a, b, lambda) |> # パラメータごとに並べ替え dplyr::mutate( dens = dgamma(x = lambda_vals, shape = a, rate = b) ) # 確率密度を計算 anime_dens_df
## # A tibble: 30,120 × 4 ## lambda a b dens ## <dbl> <dbl> <dbl> <dbl> ## 1 0 0.5 0.5 Inf ## 2 0.02 0.5 0.5 2.79 ## 3 0.04 0.5 0.5 1.96 ## 4 0.06 0.5 0.5 1.58 ## 5 0.08 0.5 0.5 1.36 ## 6 0.1 0.5 0.5 1.20 ## 7 0.12 0.5 0.5 1.08 ## 8 0.14 0.5 0.5 0.994 ## 9 0.16 0.5 0.5 0.921 ## 10 0.18 0.5 0.5 0.859 ## # … with 30,110 more rows
並べて比較する値をa_valsに指定します。
値の間隔が一定になるようにb_valsを作成します。
lambda_vals・a_vals・b_valsの全ての組み合わせを作成して、それぞれ確率密度を計算します。
パラメータごとに画面分割したアニメーション(gif画像)を作成します。
# ガンマ分布のアニメーションを作図 anime_dens_graph <- ggplot(data = anime_dens_df, mapping = aes(x = lambda, y = dens, color = as.factor(b))) + # データ geom_line(show.legend = FALSE) + # 折れ線グラフ gganimate::transition_reveal(b) + # フレーム facet_wrap(. ~ a, labeller = label_bquote(a==.(a))) + # グラフの分割 coord_cartesian(ylim = c(0, 6)) + # 軸の表示範囲 labs(title = "Gamma Distribution", subtitle = "b={frame_along}", x = expression(lambda), y = "density") # ラベル # gif画像を作成 gganimate::animate(anime_dens_graph, nframes = length(b_vals)+10, end_pause = 10, fps = 10, width = 1200, height = 600)
transition_reveal()を使ってフレームを切り替えると、過去フレームのグラフが表示され続けます。
facet_wrap()に列を指定すると、その列の値ごとにグラフを分割して描画できます。
animate()のend_pause引数を指定すると、最後のフレームでグラフが一時停止(最後のグラフを指定したフレーム数表示)します。
bの影響
続いて、$a$の値を変化させ、複数の$b$を指定して、それぞれ分布を計算します。
・作図コード(クリックで展開)
# 変化する値を指定 a_vals <- seq(from = 0.5, to = 10, by = 0.5) length(a_vals) # フレーム数 # 比較する値を指定 b_vals <- c(0.5, 1, 2, 4, 8, 12) # lambdaの値を作成 lambda_vals <- seq(from = 0, to = 5, length.out = 251) # パラメータごとにガンマ分布を計算 anime_dens_df <- tidyr::expand_grid( lambda = lambda_vals, a = a_vals, b = b_vals ) |> # 全ての組み合わせを作成 dplyr::arrange(b, a, lambda) |> # パラメータごとに並べ替え dplyr::mutate( dens = dgamma(x = lambda_vals, shape = a, rate = b) ) # 確率密度を計算 anime_dens_df
## # A tibble: 30,120 × 4 ## lambda a b dens ## <dbl> <dbl> <dbl> <dbl> ## 1 0 0.5 0.5 Inf ## 2 0.02 0.5 0.5 2.79 ## 3 0.04 0.5 0.5 1.96 ## 4 0.06 0.5 0.5 1.58 ## 5 0.08 0.5 0.5 1.36 ## 6 0.1 0.5 0.5 1.20 ## 7 0.12 0.5 0.5 1.08 ## 8 0.14 0.5 0.5 0.994 ## 9 0.16 0.5 0.5 0.921 ## 10 0.18 0.5 0.5 0.859 ## # … with 30,110 more rows
値の間隔が一定になるようにa_valsを作成します。
並べて比較する値をb_valsに指定してします。
lambda_vals・a_vals・b_valsの全ての組み合わせを作成して、それそれ確率密度を計算します。
アニメーションを作成します。
# ガンマ分布のアニメーションを作図 anime_dens_graph <- ggplot(data = anime_dens_df, mapping = aes(x = lambda, y = dens, color = as.factor(a))) + # データ geom_line(show.legend = FALSE) + # 折れ線グラフ gganimate::transition_reveal(a) + # フレーム facet_wrap(. ~ b, labeller = label_bquote(b==.(b))) + # グラフの分割 coord_cartesian(ylim = c(0, 4)) + # 軸の表示範囲 labs(title = "Gamma Distribution", subtitle = "a={frame_along}", x = expression(lambda), y = "density") # ラベル # gif画像を作成 gganimate::animate(anime_dens_graph, nframes = length(a_vals)+10, end_pause = 10, fps = 10, width = 1200, height = 600)
「aの影響」のaとbを置き換えて処理します。
パラメータと統計量の関係をアニメーションで可視化
ここまでは、パラメータと分布の関係を確認しました。次は、パラメータと統計量と歪度・尖度の関係をアニメーションで確認します。
aの影響
$a$の値を変化させ、$b$を固定して、それぞれ歪度と尖度を計算し、フレーム切替用のラベルを作成します。
・作図コード(クリックで展開)
# パラメータとして利用する値を作成 a_vals <- seq(from = 1, to = 10, by = 0.1) length(a_vals) # フレーム数 # 固定するパラメータを指定 b <- 2.5 # 作図用のlambdaの点を作成 lambda_vals <- seq(from = 0, to = 5, length.out = 501) # 歪度を計算 skewness_vec <- 2 / sqrt(a_vals) # 尖度を計算 kurtosis_vec <- 6 / a_vals # ラベル用のテキストを作成 param_vec <- paste0( "a=", a_vals, ", b=", b, ", skewness=", round(skewness_vec, 3), ", kurtosis=", round(kurtosis_vec, 3) ) head(param_vec)
## [1] "a=1, b=2.5, skewness=2, kurtosis=6" ## [2] "a=1.1, b=2.5, skewness=1.907, kurtosis=5.455" ## [3] "a=1.2, b=2.5, skewness=1.826, kurtosis=5" ## [4] "a=1.3, b=2.5, skewness=1.754, kurtosis=4.615" ## [5] "a=1.4, b=2.5, skewness=1.69, kurtosis=4.286" ## [6] "a=1.5, b=2.5, skewness=1.633, kurtosis=4"
パラメータごとに歪度と尖度を計算して、フレーム切替用のラベルとして利用します。
パラメータごとに分布を計算します。
# パラメータごとにガンマ分布を計算 anime_dens_df <- tidyr::expand_grid( lambda = lambda_vals, a = a_vals ) |> # 全ての組み合わせを作成 dplyr::arrange(a, lambda) |> # パラメータごとに並べ替え dplyr::mutate( density = dgamma(x = lambda, shape = a, rate = b), parameter = rep(param_vec, each = length(lambda_vals)) |> factor(levels = param_vec) # フレーム切替用ラベル ) # 確率密度を計算 anime_dens_df
## # A tibble: 45,591 × 4 ## lambda a density parameter ## <dbl> <dbl> <dbl> <fct> ## 1 0 1 2.5 a=1, b=2.5, skewness=2, kurtosis=6 ## 2 0.01 1 2.44 a=1, b=2.5, skewness=2, kurtosis=6 ## 3 0.02 1 2.38 a=1, b=2.5, skewness=2, kurtosis=6 ## 4 0.03 1 2.32 a=1, b=2.5, skewness=2, kurtosis=6 ## 5 0.04 1 2.26 a=1, b=2.5, skewness=2, kurtosis=6 ## 6 0.05 1 2.21 a=1, b=2.5, skewness=2, kurtosis=6 ## 7 0.06 1 2.15 a=1, b=2.5, skewness=2, kurtosis=6 ## 8 0.07 1 2.10 a=1, b=2.5, skewness=2, kurtosis=6 ## 9 0.08 1 2.05 a=1, b=2.5, skewness=2, kurtosis=6 ## 10 0.09 1 2.00 a=1, b=2.5, skewness=2, kurtosis=6 ## # … with 45,581 more rows
これまでと同様に処理します。フレーム切替用のラベル(parameter列)の値にはparam_vecを使います。
パラメータごとに統計量を計算します。
# パラメータごとに統計量を計算 anime_stat_df <- tibble::tibble( mean = a_vals / b, # 期待値 sd = sqrt(a_vals / b^2), # 標準偏差 mode = (a_vals - 1) / b, # 最頻値 parameter = factor(param_vec, levels = param_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: 364 × 3 ## parameter type statistic ## <fct> <chr> <dbl> ## 1 a=1, b=2.5, skewness=2, kurtosis=6 mean 0.4 ## 2 a=1, b=2.5, skewness=2, kurtosis=6 mode 0 ## 3 a=1, b=2.5, skewness=2, kurtosis=6 sd 0 ## 4 a=1, b=2.5, skewness=2, kurtosis=6 sd 0.8 ## 5 a=1.1, b=2.5, skewness=1.907, kurtosis=5.455 mean 0.44 ## 6 a=1.1, b=2.5, skewness=1.907, kurtosis=5.455 mode 0.0400 ## 7 a=1.1, b=2.5, skewness=1.907, kurtosis=5.455 sd 0.0205 ## 8 a=1.1, b=2.5, skewness=1.907, kurtosis=5.455 sd 0.860 ## 9 a=1.2, b=2.5, skewness=1.826, kurtosis=5 mean 0.48 ## 10 a=1.2, b=2.5, skewness=1.826, kurtosis=5 mode 0.08 ## # … with 354 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_dens_graph <- ggplot() + geom_line(data = anime_dens_df, mapping = aes(x = lambda, 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(ylim = c(0, 2.5)) + # 軸の表示範囲 labs(title = "Gamma Distribution", subtitle = "{current_frame}", x = expression(lambda), y = "density") # ラベル # gif画像を作成 gganimate::animate(anime_dens_graph, nframes = length(a_vals), fps = 10, width = 800, height = 600)
bの影響
続いて、$a$を固定し、$b$の値を変化させて、それぞれ歪度と尖度を計算し、フレーム切替用のラベルを作成します。
・作図コード(クリックで展開)
# 固定するパラメータを指定 a <- 2.5 # パラメータとして利用する値を作成 b_vals <- seq(from = 1, to = 10, by = 0.1) length(b_vals) # フレーム数 # 歪度を計算 skewness <- 2 / sqrt(a) # 尖度を計算 kurtosis <- 6 / a # ラベル用のテキストを作成 param_vec <- paste0( "a=", a, ", b=", b_vals, ", skewness=", round(skewness, 3), ", kurtosis=", round(kurtosis, 3) ) head(param_vec)
## [1] "a=2.5, b=1, skewness=1.265, kurtosis=2.4" ## [2] "a=2.5, b=1.1, skewness=1.265, kurtosis=2.4" ## [3] "a=2.5, b=1.2, skewness=1.265, kurtosis=2.4" ## [4] "a=2.5, b=1.3, skewness=1.265, kurtosis=2.4" ## [5] "a=2.5, b=1.4, skewness=1.265, kurtosis=2.4" ## [6] "a=2.5, b=1.5, skewness=1.265, kurtosis=2.4"
「aの影響」のa_valsをa、bをb_valsに置き換えると処理できます。ただし、歪度と尖度の計算に$b$は使いません。
パラメータごとに分布を計算します。
# パラメータごとにガンマ分布を計算 anime_dens_df <- tidyr::expand_grid( lambda = lambda_vals, b = b_vals ) |> # 全ての組み合わせを作成 dplyr::arrange(a, b, lambda) |> # パラメータごとに並べ替え dplyr::mutate( density = dgamma(x = lambda, shape = a, rate = b), parameter = rep(param_vec, each = length(lambda_vals)) |> factor(levels = param_vec) # フレーム切替用ラベル ) # 確率密度を計算 anime_dens_df
## # A tibble: 45,591 × 4 ## lambda b density parameter ## <dbl> <dbl> <dbl> <fct> ## 1 0 1 0 a=2.5, b=1, skewness=1.265, kurtosis=2.4 ## 2 0.01 1 0.000745 a=2.5, b=1, skewness=1.265, kurtosis=2.4 ## 3 0.02 1 0.00209 a=2.5, b=1, skewness=1.265, kurtosis=2.4 ## 4 0.03 1 0.00379 a=2.5, b=1, skewness=1.265, kurtosis=2.4 ## 5 0.04 1 0.00578 a=2.5, b=1, skewness=1.265, kurtosis=2.4 ## 6 0.05 1 0.00800 a=2.5, b=1, skewness=1.265, kurtosis=2.4 ## 7 0.06 1 0.0104 a=2.5, b=1, skewness=1.265, kurtosis=2.4 ## 8 0.07 1 0.0130 a=2.5, b=1, skewness=1.265, kurtosis=2.4 ## 9 0.08 1 0.0157 a=2.5, b=1, skewness=1.265, kurtosis=2.4 ## 10 0.09 1 0.0186 a=2.5, b=1, skewness=1.265, kurtosis=2.4 ## # … with 45,581 more rows
param_vecを使って、これまでと同様に処理します。
パラメータごとに統計量を計算します。
# パラメータごとに統計量を計算 anime_stat_df <- tibble::tibble( mean = a / b_vals, # 期待値 sd = sqrt(a / b_vals^2), # 標準偏差 mode = (a - 1) / b_vals, # 最頻値 parameter = factor(param_vec, levels = param_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: 364 × 3 ## parameter type statistic ## <fct> <chr> <dbl> ## 1 a=2.5, b=1, skewness=1.265, kurtosis=2.4 mean 2.5 ## 2 a=2.5, b=1, skewness=1.265, kurtosis=2.4 mode 1.5 ## 3 a=2.5, b=1, skewness=1.265, kurtosis=2.4 sd 0.919 ## 4 a=2.5, b=1, skewness=1.265, kurtosis=2.4 sd 4.08 ## 5 a=2.5, b=1.1, skewness=1.265, kurtosis=2.4 mean 2.27 ## 6 a=2.5, b=1.1, skewness=1.265, kurtosis=2.4 mode 1.36 ## 7 a=2.5, b=1.1, skewness=1.265, kurtosis=2.4 sd 0.835 ## 8 a=2.5, b=1.1, skewness=1.265, kurtosis=2.4 sd 3.71 ## 9 a=2.5, b=1.2, skewness=1.265, kurtosis=2.4 mean 2.08 ## 10 a=2.5, b=1.2, skewness=1.265, kurtosis=2.4 mode 1.25 ## # … with 354 more rows
「aの影響」のa_valsをa、bをb_valsに置き換えると処理できます。
「aの影響」のコードで作図できます。フレーム数はb_valsの要素数です。
# gif画像を作成 gganimate::animate(anime_dens_graph, nframes = length(a_vals), fps = 10, width = 800, height = 600)
この記事では、1次元ガウス分布の作図を確認しました。次は、乱数を生成します。
参考文献
- 須山敦志『ベイズ推論による機械学習入門』(機械学習スタートアップシリーズ)杉山将監修,講談社,2017年.
おわりに
統計量の導出もやりたいけどしばらくお預けかな。
1月28日はモーニング娘。のメジャーデビュー日です！
24周年って凄すぎる。私がハマってからは4年半くらい、これからも楽しい日々が続きますよーに。
- 2022.07.20：加筆修正しました。
【次の内容】
