はじめに

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

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

【前の内容】

www.anarchive-beta.com

【他の内容】

www.anarchive-beta.com

【今回の内容】

はじめに
負の二項分布の計算
参考文献
おわりに

負の二項分布の計算

　負の二項分布(Negative Binomial distribution)に関する計算を実装します。この記事では、Rを利用して計算します。
　負の二項分布については「負の二項分布の定義式 - からっぽのしょこ」、Pythonを利用する場合は「Python版」を参照してください。

確率の計算

　まずは、負の二項分布に従う確率(probability)を計算します。

パラメータの設定

　負の二項分布のパラメータ $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 <- 3
x

[1] 3

　失敗回数(非負の整数) $x \in \{0, 1, 2, \cdots\}$ を指定します。

　また、計算式によって必要になるパラメータ $q$ や変数 $y$ を設定します。

　成功確率から失敗確率を作成します。

# 失敗確率パラメータに変換
q <- 1 - phi
q

[1] 0.6

　成功確率 $\phi$ を用いて、失敗確率 $q = 1 - \phi$ を計算します。

　または、失敗確率から成功確率を作成します。

# 失敗確率パラメータを指定
q <- 0.4

# 失敗確率パラメータに変換
phi <- 1 - q
phi; q

[1] 0.6
[1] 0.4

　失敗確率(0から1の値) $0 \leq q \leq 1$ を指定して、成功確率 $\phi = 1 - q$ を計算します。

　失敗回数と成功回数から試行回数を作成します。

# 試行回数に変換
y <- x + r
y

[1] 8

　失敗回数 $x$ と成功回数 $r$ を用いて、試行回数 $y = x + r$ を計算します。

　または、試行回数から失敗回数を作成します。

# 確率変数の値を指定
y <- 8

# 失敗回数に変換
x <- y - r
x; y

[1] 3
[1] 8

　試行回数(正の整数) $y \in \{0, 1, 2, \cdots\}$ を指定して、失敗回数 $x = y - r$ を計算します。

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

スクラッチ実装による計算

　成功確率・失敗確率のパラメータ、試行回数の変数に関してそれぞれ実装します。

成功確率パラメータの場合

　定義式を実装して、確率を計算します。

# 定義式により確率を計算
comb <- gamma(x + r) / gamma(x + 1) / gamma(r)
prob <- comb * phi^r * (1 - phi)^x
comb; prob

[1] 35
[1] 0.1741824

　負の二項分布は、次の式で定義されます。

$\displaystyle \mathrm{NBinomial}(x \mid M, \phi) = \frac{\Gamma(x + r)}{x! \Gamma(r)} \phi^r (1 - \phi)^x$

　階乗 $x!$ の計算は、ガンマ関数 $\Gamma(x + 1) = x!$ に置き換えて行います。ガンマ関数 $\Gamma(x)$ は、gamma() で計算できます。
　または、二項係数(組み合わせ) ${}_n\mathrm{C}_m$ を、choose() で計算します。

# 二項係数を計算
comb <- choose(n = x+r-1, k = x)
comb

[1] 35

　全体の数の引数(第1引数) n に $x + r - 1$ 、選択する数の引数(第2引数) k に $x$ を指定します。

　対数をとった定義式を実装して、確率を計算します。

# 定義式により確率を計算
log_comb <- lgamma(x + r) - lgamma(x + 1) - lgamma(r)
log_prob <- log_comb + r * log(phi) + x * log(1 - phi)
prob     <- exp(log_prob)
log_prob; prob

[1] -1.747652
[1] 0.1741824

　対数をとった負の二項分布は、次の式になります。

$\displaystyle \begin{aligned} \log \mathrm{NBinomial}(x \mid r, p) &= \log \Gamma(x + r) - \log x! - \log \Gamma(r) \\ &\qquad + r \log p + x \log (1 - p) \end{aligned}$

　指数と対数の関係 $\exp(\log x) = x$ より、計算結果の指数をとると確率が得られます。

$\displaystyle \mathrm{NBinomial}(x \mid r, \phi) = \exp \Bigl( \log \mathrm{NBinomial}(x \mid r, \phi) \Bigr)$

　自然対数 $\log x$ は log()、対数ガンマ関数 $\log \Gamma(x)$ は lgamma()、指数関数 $\exp(x)$ は exp() で計算できます。
　または、対数二項係数 $\log({}_n\mathrm{C}_m)$ を、lchoose() で計算します。

# 二項係数を計算
log_comb <- lchoose(n = x+r-1, k = x)
comb     <- exp(log_comb)
log_comb; comb

[1] 3.555348
[1] 35

　choose() と同様に引数を指定します。

　ガンマ関数は、階乗の計算なので値が大きくなり、発散することがあります。

# 発散する例
gamma(170:173)

[1] 4.269068e+304 7.257416e+306           Inf           Inf

　対数ガンマ関数を使うことで、発散を回避して計算できることもあります。

# 対策例
lgamma(170:173)

[1] 701.4373 706.5731 711.7147 716.8622

失敗確率パラメータの場合

　続いて、失敗確率 $q$ を用いて、確率を計算します。

# 定義式により確率を計算
comb <- gamma(x + r) / gamma(x + 1) / gamma(r)
prob <- comb * (1 - q)^r * q^x
comb; prob

[1] 35
[1] 0.1741824

　パラメータに失敗確率 $q$ を用いると、次の式になります。

$\displaystyle \mathrm{NBinomial}(x \mid r, 1-q) = \frac{\Gamma(x + r)}{x! \Gamma(r)} (1 - q)^r q^x$

　対数をとった定義式を実装して、確率を計算します。

# 定義式により確率を計算
log_comb <- lgamma(x + r) - lgamma(x + 1) - lgamma(r)
log_prob <- log_comb + r * log(1 - q) + x * log(q)
prob     <- exp(log_prob)
log_prob; prob

[1] -1.747652
[1] 0.1741824

　対数をとると、次の式になります。

$\displaystyle \begin{aligned} \log \mathrm{NBinomial}(x \mid r, 1-q) &= \log \Gamma(x + r) - \log x! - \log \Gamma(r) \\ &\qquad + r \log (1 - q) + x \log q \end{aligned}$

試行回数変数の場合

　試行回数 $y$ を用いて、確率を計算します。

# 定義式により確率を計算
comb <- gamma(y) / gamma(y - r + 1) / gamma(r)
prob <- comb * phi^r * (1 - phi)^(y - r)
comb; prob

[1] 35
[1] 0.1741824

　確率変数に試行回数 $y$ を用いると、次の式になります。

$\displaystyle \mathrm{NBinomial}(y \mid r, p) = \frac{\Gamma(y)}{(y - r)! \Gamma(r)} p^r (1 - p)^{y-r}$

　対数をとった定義式を実装して、確率を計算します。

# 定義式により確率を計算
log_comb <- lgamma(y) - lgamma(y - r + 1) - lgamma(r)
log_prob <- log_comb + r * log(phi) + (y - r) * log(1 - phi)
prob     <- exp(log_prob)
log_prob; prob

[1] -1.747652
[1] 0.1741824

　対数をとると、次の式になります。

$\displaystyle \begin{aligned} \log \mathrm{NBinomial}(y \mid r, p) &= \log \Gamma(y) - \log (y - r)! - \log \Gamma(r) \\ &\qquad + r \log p + (y - r) \log (1 - p) \end{aligned}$

関数による計算

　成功確率・失敗確率のパラメータ、試行回数の変数に関してそれぞれ実装します。

成功確率パラメータの場合

　組み込み関数を使って、確率を計算します。

# 関数により確率を計算
prob <- dnbinom(x = x, size = r, prob = phi)
prob

[1] 0.1741824

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

　log = TRUE を指定すると、対数をとった確率を返します。

# 関数により確率を計算
log_prob <- dnbinom(x = x, size = r, prob = phi, log = TRUE)
prob     <- exp(log_prob)
log_prob; prob

[1] -1.747652
[1] 0.1741824

　計算結果の指数をとると確率が得られます。

失敗確率パラメータの場合

　続いて、失敗確率 $q$ を用いて、確率を計算します。

# 関数により確率を計算
prob <- dnbinom(x = x, size = r, prob = 1-q)
prob

[1] 0.1741824

　成功確率 $\phi = 1-q$ に変換して prob 引数に指定します。

試行回数変数の場合

　試行回数 $y$ を用いて、確率を計算します。

# 関数により確率を計算
prob <- dnbinom(x = y-r, size = r, prob = phi)
prob

[1] 0.1741824

　失敗回数 $x = y-r$ に変換して k 引数に指定します。

　ここまでは、負の二項分布の確率の計算を確認しました。次からは、分布の特徴などを表す数値の計算を確認していきます。

スポンサードリンク

統計量の計算

　次は、負の二項分布の統計量(statistics)を計算します。
　計算式については「統計量を導出」を参照してください。

　計算式を実装して、期待値を計算します。

# 期待値を計算
E_x <- r * (1 - phi) / phi
E_x

[1] 3.333333

　負の二項分布の期待値は、次の式で計算できます。

$\displaystyle \mathbb{E}[x] = \frac{r (1 - \phi)}{\phi}$

　分散を計算します。

# 分散を計算
V_x <- r * (1 - phi) / phi^2
V_x

[1] 5.555556

　分散は、次の式で計算できます。

$\displaystyle \mathbb{V}[x] = \frac{r (1 - \phi)}{\phi^2}$

　最頻値を計算します。

# 最頻値を計算
mode_x <- ifelse(
  test = r > 1, 
  yes  = floor((r - 1) * (1 - phi) / phi), 
  no   = 1
)
mode_x

[1] 2

　最頻値(モード)は、次の式で計算できます。

$\displaystyle \mathrm{mode}[x] = \begin{cases} \left\lfloor \frac{(r - 1) (1 - \phi)}{\phi} \right\rfloor &\quad (r \gt 1) \\ 0 &\quad (r \leq 1) \end{cases}$

　床関数 $\lfloor x \rfloor$ は、floor() で計算できます。

モーメントの計算

　最後に、負の二項分布の歪度(skewness)と尖度(kurtosis)を計算します。
　詳しくはいつか書きたい。

　計算式を実装して、歪度を計算します。

# 歪度を計算
skew <- (2 - phi) / sqrt((1 - phi) * r)
skew

[1] 0.9899495

　負の二項分布の歪度は、次の式で計算できます。

$\displaystyle \mathrm{Skewness} = \frac{ 2 - \phi }{ \sqrt{(1 - \phi) r} }$

　ルート $\sqrt{x}$ は、sqrt() で計算できます。

　尖度を計算します。

# 尖度を計算
kurt <- 6 / r + phi^2 / (1 - phi) / r
kurt

[1] 1.38

　尖度は、次の式で計算できます。

$\displaystyle \mathrm{Kurtosis} = \frac{6}{r} + \frac{ \phi^2 }{ (1 - \phi) r }$

　この記事では、負の二項分布に関する計算を確認しました。次の記事では、確率分布のグラフを作成します。

参考文献

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

作者:須山敦志
講談社

Amazon

おわりに

　他のパラメータとの兼ね合いから失敗確率の表記を $\psi = 1 - \phi$ にしたかったのですが、一度切りの記号でギリシャ文字が登場しても思考の邪魔だと判断して、 $q$ にしました。次点で成功確率を $p$ にするという手も考えましたが、他の分布との兼ね合いから止めておきました。

　2025年10月19日は、BEYOOOOONDSの結成7周年の日です。

　最KOOOOO！！！！！！！

【次の内容】

からっぽのしょこ

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

【R】負の二項分布の計算

はじめに

負の二項分布の計算

確率の計算

パラメータの設定

スクラッチ実装による計算

成功確率パラメータの場合

失敗確率パラメータの場合

試行回数変数の場合

関数による計算

成功確率パラメータの場合

失敗確率パラメータの場合

試行回数変数の場合

統計量の計算

モーメントの計算

参考文献

おわりに