はじめに

　機械学習や統計学で登場する各種の確率分布について、「計算式の導出・計算のスクラッチ実装・計算過程や結果の可視化」などの「数式・プログラム・図」を用いた解説により、様々な角度から理解を目指すシリーズです。

　この記事では、負の二項分布のグラフ作成ついてR言語を使って確認します。

【前の内容】

www.anarchive-beta.com

【他の内容】

www.anarchive-beta.com

【今回の内容】

はじめに
負の二項分布の作図
- グラフの作成
- 複数グラフの作成
  - パラメータの設定
  - 確率分布の作図
    - 重ねて描画
    - 分割して描画
参考文献
おわりに

負の二項分布の作図

　負の二項分布(Negative Binomial distribution)のグラフを作成します。この記事では、Rのggplot2パッケージを利用して作図します。
　負の二項分布については「負の二項分布の定義式 - からっぽのしょこ」、確率計算については「【R】負の二項分布の計算 - からっぽのしょこ」、Pythonを利用する場合は「Python版」を参照してください。

　利用するパッケージを読み込みます。

# 利用パッケージ
library(tidyverse)

　この記事では、基本的に パッケージ名::関数名() の記法を使うので、パッケージを読み込む必要はありません。ただし、作図コードについては(ごちゃごちゃしないように)パッケージ名を省略するため、ggplot2 を読み込む必要があります。
　また、ネイティブパイプ(baseパイプ)演算子 |> を使います。(基本的には)magrittrパッケージのパイプ演算子 %>% に置き換えても処理できますが、その場合は magrittr を読み込む必要があります。

グラフの作成

　まずは、負の二項分布のグラフを作成します。

パラメータの設定

　負の二項分布のパラメータ $r, \phi$ を設定します。

# 成功回数パラメータを指定
r <- 5

# 成功確率パラメータを指定
phi <- 0.4
r; phi

[1] 5
[1] 0.4

　成功回数(正の整数) $r \in \{1, 2, \cdots\}$ 、成功確率(0から1の値) $0 \leq \phi \leq 1$ を指定します。

　確率変数 $x$ の作図範囲を設定します。

# x軸の範囲を設定
x_min <-0
#x_max <- 25
u <- 5
x_max <- (r * (1-phi) / phi) |> # 基準値を指定
  (\(.) {. * 3})() |> # 倍率を指定
  (\(.) {ceiling(. /u)*u})() # u単位で切り上げ

# x軸の値を作成
x_vec <- seq(from = x_min, to = x_max, by = 1)
x_min; x_max; head(x_vec)

[1] 0
[1] 25
[1] 0 1 2 3 4 5

　失敗回数(非負の整数) $x \in \{0, 1, 2, \ldots\}$ の数値ベクトルを作成します。
　この例では、指定したパラメータ(から求める変数の期待値 $\mathbb{E}[x] = \frac{r (1 - \phi)}{\phi}$ ) r, phi を使って、範囲を設定しています。

　続いて、設定した値に従う確率を計算して、グラフを作成していきます。

確率分布の計算

　 $x$ の値ごとに負の二項分布に従う確率 $\mathrm{NB}(x \mid r, \phi)$ を計算します。

# 確率を計算
prob_df <- tidyr::tibble(
  x    = x_vec, # 確率変数
  prob = dnbinom(x = x, size = r, prob = phi) # 確率
)
prob_df

# A tibble: 26 × 2
       x   prob
   <dbl>  <dbl>
 1     0 0.0102
 2     1 0.0307
 3     2 0.0553
 4     3 0.0774
 5     4 0.0929
 6     5 0.100 
 7     6 0.100 
 8     7 0.0946
 9     8 0.0851
10     9 0.0738
# ℹ 16 more rows

　負の二項分布の確率は、dnbinom() で計算できます。確率変数の引数 x に $x$ (のベクトルや列)、成功回数の引数 size に $r$ 、成功確率の引数 prob に $\phi$ を指定します。

　確率変数 $x$ を x 列、 $x$ の確率 $p(x \mid r, \phi)$ を prob 列として、データフレームに格納します。

確率分布の作図

　パラメータラベルを作成します。装飾用の処理です。

# ラベル用の文字列を作成
param_lbl <- paste0(
  "list(", 
  "r == ", r, ", ", 
  "phi == ", round(phi, digits = 2), 
  ")"
) |> 
  parse(text = _)
param_lbl

expression(list(r == 5, phi == 0.4))

　ギリシャ文字などの記号を使った数式を表示する場合は、expression() の記法を使います。等号は "=="、複数の(数式上の)変数は list(変数1, 変数2) で並べて描画できます。
　オブジェクト(プログラム上の変数)の値を使う場合は、文字列を作成して parse() の text 引数に渡します。

　負の二項分布の棒グラフを作成します。

# 負の二項分布を作図:棒グラフ
ggplot() + 
  geom_bar(
    data = prob_df, 
    mapping = aes(x = x, y = prob), 
    stat = "identity", #position = "identity", 
    fill = "#00A968"
  ) + # バー
  scale_x_continuous(breaks = x_vec, minor_breaks = FALSE) + # x軸目盛
  labs(
    title = "Negative Binomial distribution", 
    subtitle = param_lbl, 
    x = "x", y = "probability"
  )

　棒グラフは、geom_bar() で描画できます。縦軸の値を直接指定する(データフレームに計算済み値が格納されている)場合は、stat 引数に "identity" を指定します。

　幹図を作成します。

# 負の二項分布を作図:幹図
ggplot() + 
  geom_segment(
    data = prob_df, 
    mapping = aes(x = x, y = 0, xend = x, yend = prob), 
    linewidth = 1
  ) + # 線分
  geom_point(
    data = prob_df, 
    mapping = aes(x = x, y = prob), 
    color = "#00A968", size = 5
  ) + # 点
  scale_x_continuous(breaks = x_vec, minor_breaks = FALSE) + # x軸目盛
  labs(
    title = "Negative Binomial distribution", 
    subtitle = param_lbl, 
    x = "x", y = "probability"
  )

　幹図は、geom_segment() と geom_point() を使って描画できます。geom_segment() の x, xend 引数にx座標、y, yend 引数に0と確率値を指定して、線分を描画します。先端に重ねるように、geom_point() で確率値の位置に点を描画します。

　折れ線グラフを作成します。

# 負の二項分布を作図:折れ線グラフ
ggplot() + 
  geom_line(
    data = prob_df, 
    mapping = aes(x = x, y = prob), 
    linewidth = 1
  ) + # 折れ線
  geom_point(
    data = prob_df, 
    mapping = aes(x = x, y = prob), 
    color = "#00A968", size = 5
  ) + # 点
  scale_x_continuous(breaks = x_vec, minor_breaks = FALSE) + # x軸目盛
  labs(
    title = "Negative Binomial distribution", 
    subtitle = param_lbl, 
    x = "x", y = "probability"
  )

　折れ線グラフは、geom_line() で描画できます。
　負の二項分布は離散型の確率分布なので、変数がとり得る値(離散値)の位置に点を描画します。

　以上で、基本的な作図処理を確認しました。

スポンサードリンク

複数グラフの作成

　次は、複数のパラメータの負の二項分布を並べて描画します。

パラメータの設定

　複数個のパラメータ $r, \phi$ を設定します。

# 成功回数パラメータを指定
r_vals <- c(1, 5, 10)

# 成功確率パラメータを指定
phi_vals <- c(0.5, 0.4, 0.3)
r_vals; phi_vals

[1]  1  5 10
[1] 0.5 0.4 0.3

　成功回数 $r_i$ 、成功確率 $\phi_i$ をそれぞれ数値ベクトルで指定します。

確率分布の作図

　1つの領域に複数の分布を重ねて表示する図と、複数の領域に分割してそれぞれの分布を表示する図の2パターンを作成します。

重ねて描画

　こちらのの図は、r_vals と phi_vals が同じパラメータ数である必要があります。

　確率変数 $x$ の作図範囲を設定します。

# x軸の範囲を指定
x_min <- 0
#x_max <- 50
u <- 5
x_max <- (r_vals * (1-phi_vals) / phi_vals) |> # 基準値を指定
  max() |> 
  (\(.) {. * 2})() |> # 倍率を指定
  (\(.) {ceiling(. /u)*u})() # u単位で切り上げ

# x軸の値を作成
x_vec <- seq(from = x_min, to = x_max, by = 1)
x_min; x_max; head(x_vec)

[1] 0
[1] 50
[1] 0 1 2 3 4 5

　「グラフの作成」のときと同様に処理します。
　この例では、指定したパラメータ(から求める変数の期待値 $\mathbb{E}[x] = \frac{r_i (1 - \phi_i)}{\phi_i}$ の最大値) r_vals, lambda_vals を使って、範囲を設定しています。

　 $r, \phi$ と $x$ の組み合わせごとに確率を計算します。

# パラメータごとに確率を計算
res_prob_df <- tibble::tibble(
  r   = r_vals,  # 成功回数
  phi = phi_vals # 成功確率
) |> 
  dplyr::mutate(
    i = dplyr::row_number() # パラメータ番号
  ) |> 
  tidyr::expand_grid(
    x = x_vec # 確率変数
  ) |> # パラメータごとに変数を複製
  dplyr::mutate(
    prob = dnbinom(x = x, size = r, prob = phi) # 確率
  ) |> 
  dplyr::select(i, r, phi, x, prob)
res_prob_df

# A tibble: 153 × 5
       i     r   phi     x     prob
   <int> <dbl> <dbl> <dbl>    <dbl>
 1     1     1   0.5     0 0.5     
 2     1     1   0.5     1 0.25    
 3     1     1   0.5     2 0.125   
 4     1     1   0.5     3 0.0625  
 5     1     1   0.5     4 0.0312  
 6     1     1   0.5     5 0.0156  
 7     1     1   0.5     6 0.00781 
 8     1     1   0.5     7 0.00391 
 9     1     1   0.5     8 0.00195 
10     1     1   0.5     9 0.000977
# ℹ 143 more rows

　パラメータ $r_i, \phi_i$ を格納して、変数 $x$ との全ての組み合わせを expand_grid() で作成します。パラメータのペアと変数の組み合わせ(行)ごとに dnbinom() で負の二項分布の確率を計算します。
　パラメータのペア(分布)ごとに色分けする場合は、パラメータ番号(通し番号)を row_number() で割り振っておきます。

　パラメータ番号 $i$ を i 列、成功回数 $r_i$ を r 列、成功確率 $\phi_i$ を phi 列、確率変数 $x$ を x 列、 $x$ の確率 $p(x \mid r_i, \phi_i)$ を prob 列として、データフレームに格納します。

　パラメータラベルを作成します。

# ラベル用の文字列を作成
param_lbl_vec <- paste0(
  "list(", 
  "r == ", r_vals, ", ", 
  "phi == ", round(phi_vals, digits = 2), 
  ")"
) |> 
  parse(text = _)
param_lbl_vec[1:3]

expression(list(r == 1, phi == 0.5), list(r == 5, phi == 0.4), 
    list(r == 10, phi == 0.3))

　「グラフの作成」のときと同様に処理して、ラベル用のベクトルを作成します。

　パラメータごとの負の二項分布を重ねて描画します。

# パラメータごとに負の二項分布を作図
ggplot() + 
  geom_bar(
    data = res_prob_df, 
    mapping = aes(x = x, y = prob, fill = factor(i)), 
    stat = "identity", position = "identity", 
    alpha = 0.5
  ) + # バー
  geom_line(
    data = res_prob_df, 
    mapping = aes(x = x, y = prob, color = factor(i)), 
    linewidth = 1
  ) + # 折れ線
  geom_point(
    data = res_prob_df, 
    mapping = aes(x = x, y = prob, color = factor(i)), 
    size = 3
  ) + # 点
  #scale_x_continuous(breaks = x_vec, minor_breaks = FALSE) + # x軸目盛
  scale_fill_hue(labels = param_lbl_vec) + # (数式表示用)
  scale_color_hue(labels = param_lbl_vec) + # (数式表示用)
  labs(
    title = "Negative Binomial distribution", 
    fill = "parameter", color = "parameter", 
    x = expression(x), 
    y = "probability"
  )

　「グラフの作成」のときと同様に各種グラフを描画します。
　複数カテゴリの棒グラフを描画する場合は、position 引数に "identity" を指定します。
　パラメータ(任意のカテゴリ)ごとに色付けする場合は、因子型のパラメータ番号(カテゴリを示すIDなどの列)を fill, color 引数に指定します。数値型の列を指定すると値の大小に応じて配色されます。線の種類などで描き分ける場合も同様に指定します。
　凡例に数式を表示する場合は、scale_***_hue() などの labels 引数に expression() 記法で指定します。

分割して描画

　こちらの図は、r_vals と phi_vals が異なるパラメータ数でも処理できます。

　確率変数 $x$ の作図範囲を設定します。

# x軸の範囲を指定
x_min <- 0
#x_max <- 50
u <- 5
x_max <- outer(
  X   = r_vals, 
  Y   = phi_vals, 
  FUN = \(r, p) {r * (1-p) / p}
) |> # 基準値を指定
  max() |> 
  (\(.) {. * 2})() |> # 倍率を指定
  (\(.) {ceiling(. /u)*u})() # u単位で切り上げ

# x軸の値を作成
x_vec <- seq(from = x_min, to = x_max, by = 1)
x_min; x_max; head(x_vec)

[1] 0
[1] 50
[1] 0 1 2 3 4 5

　「グラフの作成」のときと同様に処理します。
　この例では、指定したパラメータ(から求める変数の期待値 $\mathbb{E}[x] = \frac{r_i (1 - \phi_j)}{\phi_j}$ の最大値) r_vals, lambda_vals を使って、範囲を設定しています。

　 $r$ と $\phi$ の組み合わせごとに確率を計算します。

# パラメータごとに確率を計算
res_prob_df <- tidyr::expand_grid(
  r   = r_vals,  # 成功回数
  phi = phi_vals # 成功確率
) |> # 全ての組み合わせを作成
  dplyr::mutate(
    i = dplyr::row_number() # パラメータ番号
  ) |> 
  tidyr::expand_grid(
    x = x_vec # 確率変数
  ) |> # パラメータごとに変数を複製
  dplyr::mutate(
    prob = dnbinom(x = x, size = r, prob = phi) # 確率
  ) |> 
  dplyr::select(i, r, phi, x, prob)
res_prob_df

# A tibble: 459 × 5
       i     r   phi     x     prob
   <int> <dbl> <dbl> <dbl>    <dbl>
 1     1     1   0.5     0 0.5     
 2     1     1   0.5     1 0.25    
 3     1     1   0.5     2 0.125   
 4     1     1   0.5     3 0.0625  
 5     1     1   0.5     4 0.0312  
 6     1     1   0.5     5 0.0156  
 7     1     1   0.5     6 0.00781 
 8     1     1   0.5     7 0.00391 
 9     1     1   0.5     8 0.00195 
10     1     1   0.5     9 0.000977
# ℹ 449 more rows

　2つのパラメータ $r_i, \phi_j$ の全ての組み合わせを expand_grid() で作成して、さらに変数 $x$ との全ての組み合わせを作成して、確率を計算します。

　パラメータの組み合わせごとに描画領域を分割して分布を描画します。

# パラメータごとに負の二項分布を作図
ggplot() + 
  geom_bar(
    data = res_prob_df, 
    mapping = aes(x = x, y = prob), 
    stat = "identity", position = "identity"
  ) + # バー
  #scale_x_continuous(breaks = x_vec, minor_breaks = FALSE) + # x軸目盛
  facet_grid(
    rows = vars(r), 
    cols = vars(phi), 
    labeller = label_bquote(
      r == .(r), 
      phi == .(round(phi, digits = 2))
    )
  ) + # 描画領域の分割
  labs(
    title = "Negative Binomial distribution", 
    x = expression(x), 
    y = "probability"
  )

　facet_grid() で描画領域を分割してパラメータごとにグラフを描画します。rows 引数に縦に並べる列、cols 引数に横に並べる列を指定します。列の値(や文字列)が昇順に並べ替えられて描画されます。任意の順番に描画する場合は、列の値を因子型に変換して、因子レベルを設定しておく必要があります。
　labeller 引数で個別にラベルを指定できます。数式を表示する場合は、label_bquote() を使います。データフレームの値(要素)を使う場合は .(列名) の形式で指定します。

　この例では、 $r$ の変化を縦方向、 $\phi$ の変化を横方向に並べています。

　以上で、複数の確率分布の作図処理を確認しました。

　この記事では、負の二項分布のグラフを作成しました。次の記事では、パラメータと分布の関係を可視化します。