はじめに

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

　この記事では、負の二項分布の乱数生成についてR言語を使って確認します。

【前の内容】

www.anarchive-beta.com

【他の内容】

www.anarchive-beta.com

【今回の内容】

はじめに
負の二項分布の乱数生成
- サンプリング
- サンプルサイズの影響
  - サンプルサイズと形状の関係
参考文献
おわりに

負の二項分布の乱数生成

　負の二項分布(Negative Binomial distribution)に従う乱数(random number)を生成して、グラフやアニメーションで確認します。この記事では、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$ を指定します。

　続いて、設定した値に従う乱数を生成して、グラフを作成していきます。

乱数の生成

　サンプルサイズ $N$ を指定して、負の二項分布の乱数 $x_n$ を生成します。

# サンプルサイズを指定
N <- 100

# 負の二項分布の乱数を生成
x_n <- rnbinom(n = N, size = r, prob = phi)
head(x_n)

[1] 5 3 4 8 7 8

　サンプルサイズ(データ数) $N$ を指定します。
　 $N$ 個のサンプル(乱数)の集合 $\mathbf{x} = \{x_1, \cdots, x_N\}$ を作成します。各要素は、失敗回数(非負の整数) $x_n \in \{0, 1, 2, \cdots\}$ です。
　負の二項分布の乱数は、rnbinom() で生成できます。サンプルサイズの引数 n に $N$ 、成功回数の引数 size に $r$ 、成功確率の引数 prob に $\phi$ を指定します。

　以上で、乱数を作成できました。続いて、乱数を可視化していきます。

乱数の集計

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

# x軸の範囲を設定
u <- 5
x_max <- x_n |> # 基準値を指定
  max() |> 
  (\(.) {ceiling(. /u)*u})() # u単位で切り上げ
#x_max <- 25

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

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

　集計や作図の処理用に、単位時間における失敗回数(非負の整数) $x \in \{0, 1, 2, \cdots\}$ の数値ベクトルを作成します。
　この例では、生成したサンプル x_n の最大値を使って、範囲を設定しています。

　サンプル集合 $\mathbf{x}$ を集計します。

# 乱数を集計
freq_df <- tibble::tibble(
  x = x_n # サンプル値
) |> 
  dplyr::count(
    x, name = "freq" # 度数
  ) |> 
  dplyr::mutate(
    rel_freq = freq / N # 相対度数
  )
freq_df

# A tibble: 19 × 3
       x  freq rel_freq
   <int> <int>    <dbl>
 1     1     3     0.03
 2     2     8     0.08
 3     3     8     0.08
 4     4    13     0.13
 5     5    11     0.11
 6     6     6     0.06
 7     7     8     0.08
 8     8     8     0.08
 9     9     5     0.05
10    10     7     0.07
11    11     5     0.05
12    12     6     0.06
13    13     3     0.03
14    14     2     0.02
15    16     3     0.03
16    17     1     0.01
17    19     1     0.01
18    22     1     0.01
19    24     1     0.01

　サンプル値 $\mathbf{x}$ をデータフレームに格納して、count() で重複する値をカウントし(度数を求め)ます。
　度数を $N_x$ ( 値が $1$ の度数を $N_1$ で表す)として、サンプルサイズ $N$ で割ると相対度数 $\frac{N_x}{N}$ が求まります。
　観測されなかった変数の値もデータフレームに含める場合は、次の処理を追加します。

# 未観測値を補完
freq_df <- freq_df |> 
  tidyr::complete(
    x = x_vec, 
    fill = list(freq = 0, rel_freq = 0)
  )
freq_df

# A tibble: 26 × 3
       x  freq rel_freq
   <dbl> <int>    <dbl>
 1     0     0     0   
 2     1     3     0.03
 3     2     8     0.08
 4     3     8     0.08
 5     4    13     0.13
 6     5    11     0.11
 7     6     6     0.06
 8     7     8     0.08
 9     8     8     0.08
10     9     5     0.05
# ℹ 16 more rows

　 $x$ がとり得る値を指定して、complete() で値を補完します。欠損値となる他の列の値を fill 引数にリスト形式で指定します。

　確率変数 $x$ を x 列、 $x$ の度数 $N_x$ を freq 列、相対度数 $\frac{N_x}{N}$ を rel_freq 列として、データフレームに格納します。

乱数の作図

　負の二項分布の乱数のヒストグラムを作成します。

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

# サンプルの度数を作図
ggplot() + 
  geom_bar(
    data = freq_df, 
    mapping = aes(x = x, y = freq), 
    stat = "identity", position = "identity", width = 1, 
    fill = "#00A968"
  ) + # 度数
  scale_x_continuous(breaks = x_vec, minor_breaks = FALSE) + # x軸目盛
  labs(
    title = "Negative Binomial distribution", 
    subtitle = param_lbl, 
    x = expression(x), 
    y = "frequency"
  )

　棒グラフの描画関数 geom_bar() でヒストグラムを作成します。横軸の引数 x に階級値、縦軸の引数 y に度数を指定します。また、(棒グラフではなく)ヒストグラムの体裁にするには、バーの幅の引数 width に 1 (階級幅)を指定します。

　または、集計前のサンプルを用いて、ヒストグラムを作成します。

# サンプルを格納
sample_df <- tibble::tibble(
  x = x_n
)

# 階級幅を設定
bin_size <- 1

# 境界値を作成
center_vec <- seq(from = 0, to = x_max, by = bin_size) # 階級値
bin_vec    <- c(center_vec, max(center_vec)+bin_size) - 0.5*bin_size # 境界値

# サンプルの度数を作図
ggplot() + 
  geom_histogram(
    data = sample_df, 
    mapping = aes(x = x, y = after_stat(count)), 
    stat = "bin", position = "identity", breaks = bin_vec,
    fill = "#00A968"
  ) + # 度数
  scale_x_continuous(breaks = x_vec, minor_breaks = FALSE) + # x軸目盛
  coord_cartesian(xlim = c(-0.5*bin_size, x_max+0.5*bin_size)) + # 表示範囲
  labs(
    title = "Negative Binomial distribution", 
    subtitle = param_lbl, 
    x = expression(x), 
    y = "frequency"
  )

　ヒストグラムは、geom_histogram() で描画できます。横軸の引数 x にサンプルデータ、縦軸の引数 y に after_stat(count) を指定します。また、バーの境界の引数 breaks に境界値(階級の下限値・上限値)を指定します。

　負の二項分布に従う確率を計算します。装飾用の処理です。

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

　詳しくは「分布の作図」を参照してください。

　相対度数のグラフを作成します。

# サンプルの相対度数を作図
ggplot() + 
  geom_bar(
    data = freq_df, 
    mapping = aes(x = x, y = rel_freq, linetype = "sample"), 
    stat = "identity", position = "identity", 
    fill = "#00A968", alpha = 0.5
  ) + # 相対度数
  geom_bar(
    data = prob_df, 
    mapping = aes(x = x, y = prob, linetype = "generator"), 
    stat = "identity", position = "identity", 
    fill = NA, color = "darkgreen", linewidth = 1
  ) + # 生成確率
  #scale_x_continuous(breaks = x_vec, minor_breaks = FALSE) + # x軸目盛
  scale_y_continuous(
    sec.axis = sec_axis(transform = ~ .*N, name = "frequency")
  ) + # 度数軸目盛
  scale_linetype_manual(
    breaks = c("generator", "sample"), 
    values = c("dashed", "blank"), 
    labels = c("generator", "random number"), 
    name   = "distribution"
  ) + # 凡例の表示用
  guides(
    linetype = guide_legend(override.aes = list(linewidth = 0.5))
  ) + # 凡例の体裁
  labs(
    title = "Negative Binomial distribution", 
    subtitle = param_lbl, 
    x = expression(x), 
    y = "relative frequency, probability"
  )

　形状の比較用に、負の二項分布に従う確率(サンプルの生成分布)を破線で示します。

　y 引数に相対度数(サンプル全体に対する割合)を指定します。

　以上で、基本的な乱数生成の処理を確認しました。

スポンサードリンク

サンプルサイズの影響

　次は、サンプルサイズを増やしてグラフの変化をアニメーションで確認します。
　作図コードについては「Probability-Distribution/code/negative_binomial/random_number.R at main · anemptyarchive/Probability-Distribution · GitHub」を参照してください。

サンプルサイズと形状の関係

　サンプルサイズ $N$ を大きくしたときの負の二項分布の乱数のヒストグラム(度数)の変化をアニメーションにします。

　上の図は1サンプルずつ、下の図は複数サンプルずつ増やしたときの様子です。

　相対度数の変化をアニメーションにします。

　サンプルが増えるほど、生成分布の形状に近付くのが分かります。

　以上で、サンプルサイズの影響を確認しました。

　この記事では、負の二項分布の乱数を生成しました。

参考文献

機械学習スタートアップシリーズ　ベイズ推論による機械学習入門

作者:須山敦志
講談社

Amazon

おわりに

　ここまでの作業でPython版と合わせて離散分布・連続分布・混合分布の解説や処理の表記揺れなどを統一できたので、以降の分布の作業が捗るはず(n回目の希望的観測)。これで負の二項分布について新規でまとめられたので、(ポアソン分布に関しては)ベイズ推論の方の作業も進められるはずと思ったけど、ガンマ分布の修正作業がまだだった。

【次の内容】

からっぽのしょこ

読んだら書く！書いたら読む！同じ事は二度調べ(たく)ない

【R】負の二項分布の乱数生成