はじめに

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

　この記事では、混合ポアソン分布の確率の計算についてR言語を使って確認します。

【前の内容】

www.anarchive-beta.com

【他の内容】

www.anarchive-beta.com

【今回の内容】

はじめに
混合ポアソン分布の計算
- 確率の計算
参考文献
おわりに

混合ポアソン分布の計算

　混合ポアソン分布(mixture Poisson distribution)の計算を実装します。この記事では、Rを利用して計算します。
　ポアソン分布については「【R】ポアソン分布の計算 - からっぽのしょこ」、混合ポアソン分布については「混合ポアソン分布の定義式の導出 - からっぽのしょこ」、Pythonを利用する場合は「Python版」を参照してください。

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

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

　この記事では、基本的に組み込み関数を使うので、パッケージを読み込む必要はありません。ただし、データフレーム上で計算を行う場合には、tidyverseパッケージを利用します。その際も、パッケージ名::関数名() の記法を使うので、パッケージを読み込む必要はありません。
　また、ネイティブパイプ(baseパイプ)演算子 |> を使います。(基本的には)magrittrパッケージのパイプ演算子 %>% に置き換えても処理できますが、その場合は magrittr を読み込む必要があります。

確率の計算

　まずは、混合ポアソン分布に従う確率を計算します。

パラメータの設定

　混合ポアソン分布のパラメータ $\boldsymbol{\lambda}$ を設定します。

# クラスタごとのパラメータを指定
lambda_k <- c(5.5, 10, 16.8)
lambda_k

[1]  5.5 10.0 16.8

　各クラスタ(クラス)に対応する $K$ 個のポアソン分布のパラメータ $\boldsymbol{\lambda} = \{\lambda_1, \cdots, \lambda_K\}$ を数値ベクトルで指定します。各値は、単位時間における発生回数の期待値であり、正の値 $\lambda_k \gt 0$ を満たす必要があります。

　クラスタの混合比率 $\boldsymbol{\pi}$ を設定する。

# 混合比率を指定
pi_k <- c(0.35, 0.25, 0.4)
pi_k; sum(pi_k)

[1] 0.35 0.25 0.40
[1] 1

　カテゴリ分布のパラメータ $\boldsymbol{\pi} = (\pi_1, \cdots, \pi_K)$ を数値ベクトルで指定します。各値は、クラスタの割り当て確率であり、0から1の値 $0 \leq \pi_k \leq 1$ で、総和が1になる値 $\sum_{k=1}^K \pi_k = 1$ を満たす必要があります。

　クラスタ数 $K$ を設定します。

# クラスタ数を取得
K <- length(lambda_k)
K

[1] 3

　クラスタ数 $K$ を作成します。
　この例では、パラメータの要素数を使います。

　確率変数の実現値 $x$ を設定します。

# 確率変数の値を指定
x <- 10
x

[1] 10

　単位時間における発生回数(0以上の整数) $x \in \{0, 1, 2, \cdots\}$ をスカラで指定します。

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

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

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

## 定義式により確率を計算

# クラスタごとのポアソン分布の確率を計算
cluster_prob_k <- exp(-lambda_k) * lambda_k^x / gamma(x + 1)

# クラスタごとの重み付け確率を計算
weighted_prob_k <- pi_k * cluster_prob_k

# 周辺確率を計算
prob <- sum(weighted_prob_k)
round(cluster_prob_k, 3); round(weighted_prob_k, 3); round(prob, 3)

[1] 0.029 0.125 0.025
[1] 0.010 0.031 0.010
[1] 0.051

　混合ポアソン分布の周辺確率(全てのクラスタを考慮した生成確率)は、次の式で計算できます。

$\displaystyle \begin{aligned} p(x \mid \boldsymbol{\lambda}, \boldsymbol{\pi}) &= \sum_{k=1}^K p(x, s = k \mid \lambda_k, \boldsymbol{\pi}) \\ &= \sum_{k=1}^K p(s = k \mid \boldsymbol{\pi}) p(x \mid s = k, \lambda_k) \\ &= \sum_{k=1}^K \pi_k \mathrm{Poisson}(x \mid \lambda_k) \\ &= \sum_{k=1}^K \pi_k \exp(-\lambda_k) \frac{\lambda_k^x}{x!} \end{aligned}$

　階乗 $x!$ の計算は、ガンマ関数 $\Gamma(x + 1) = x!$ に置き換えて行います。ガンマ関数は、gamma() で計算できます。
　ネイピア数 $e$ の指数関数 $e^x = \exp(x)$ は、exp() で計算できます。

　この例では、 $K$ 個のクラスタの計算をベクトル演算で処理しています。
　クラスタ $k$ に関する総和 $\sum_{k=1}^K$ の計算をforループで処理します。

## 定義式により確率を計算

# 周辺確率を初期化
prob <- 0

# クラスタごとに処理
for(k in 1:K) {
  
  # クラスタkのパラメータを取得
  tmp_lambda <- lambda_k[k]
  tmp_pi     <- pi_k[k]
  
  # クラスタkのポアソン分布の確率を計算
  tmp_cluster_prob <- exp(-tmp_lambda) * tmp_lambda^x / gamma(x + 1)
  
  # クラスタkの重み付け確率を計算
  tmp_weighted_prob <- tmp_pi * tmp_cluster_prob
  
  # 周辺確率を計算
  prob <- prob + tmp_weighted_prob
}
round(prob, 3)

[1] 0.051

　パラメータのベクトルから順番に値(要素)を取り出して、重み付け確率を計算し prob に足していきます。始めに prob の値を 0 に初期化(値が0のオブジェクトを作成)しておく必要があります。

関数による確率の計算

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

## 関数により確率を計算

# クラスタごとのポアソン分布の確率を計算
cluster_prob_k <- dpois(x = x, lambda = lambda_k)

# クラスタごとの重み付け確率を計算
weighted_prob_k <- pi_k * cluster_prob_k

# 周辺確率を計算
prob <- sum(weighted_prob_k)
round(cluster_prob_k, 3); round(weighted_prob_k, 3); round(prob, 3)

[1] 0.029 0.125 0.025
[1] 0.010 0.031 0.010
[1] 0.051

　ポアソン分布の確率は、dpois() で計算できます。確率変数の引数 x に $x$ 、パラメータの引数 lambda に $\lambda_k$ を指定します。
　x 引数にスカラ、lambda 引数にベクトルを指定することで、複数のパラメータの計算を一度に処理できます。

　forループで計算します。

## 定義式により確率を計算

# 周辺確率を初期化
prob <- 0

# クラスタごとに処理
for(k in 1:K) {
  
  # クラスタkのパラメータを取得
  tmp_lambda <- lambda_k[k]
  tmp_pi     <- pi_k[k]
  
  # クラスタkのポアソン分布の確率を計算
  tmp_cluster_prob <- dpois(x = x, lambda = tmp_lambda)
  
  # クラスタkの重み付け確率を計算
  tmp_weighted_prob <- tmp_pi * tmp_cluster_prob
  
  # 周辺確率を計算
  prob <- prob + tmp_weighted_prob
}
round(prob, 3)

[1] 0.051

　「スクラッチ実装による確率の計算」のクラスタごとの確率の計算を関数 dpois() に置き換えて処理します。

複数データの場合

　確率変数の実現値 $\mathbf{x}$ を設定します。

# 確率変数の値を指定
x_n <- c(10, 6, 17, 13, 15)

# データ数を取得
N <- length(x_n)
N; x_n

[1] 5
[1] 10  6 17 13 15

　 $N$ 個の変数 $\mathbf{x} = \{x_1, \cdots, x_N\}$ を数値ベクトルで指定して、データ数 $N$ を作成します。

　途中計算用に、変数 $\mathbf{x}$ とパラメータ $\boldsymbol{\lambda}$ の値を複製します。

# 確率変数を複製
x_nk <- x_n |> 
  rep(times = K) |> 
  matrix(nrow = N, ncol = K)

# パラメータを複製
lambda_nk <- lambda_k |> 
  rep(each = N) |> 
  matrix(nrow = N, ncol = K)
dim(x_nk); x_nk; dim(lambda_nk); lambda_nk

[1] 5 3
     [,1] [,2] [,3]
[1,]   10   10   10
[2,]    6    6    6
[3,]   17   17   17
[4,]   13   13   13
[5,]   15   15   15
[1] 5 3
     [,1] [,2] [,3]
[1,]  5.5   10 16.8
[2,]  5.5   10 16.8
[3,]  5.5   10 16.8
[4,]  5.5   10 16.8
[5,]  5.5   10 16.8

　 $N$ 個の変数 $n = 1, \dots, N$ が行、 $K$ 個のクラスタ $k = 1, \dots, K$ が列に対応するように、 $\mathbf{x}, \boldsymbol{\lambda}$ それぞれのベクトルを複製します。

　混合ポアソン分布の確率を計算します。

# クラスタごとのポアソン分布の確率を計算
cluster_prob_nk <- dpois(x = x_nk, lambda = lambda_nk)

# クラスタごとの重み付け確率を計算
weighted_prob_nk <- t(pi_k * t(cluster_prob_nk))

# 周辺確率を計算
prob_n <- rowSums(weighted_prob_nk)
round(cluster_prob_nk, 3); round(weighted_prob_nk, 3); round(prob_n, 3)

      [,1]  [,2]  [,3]
[1,] 0.029 0.125 0.025
[2,] 0.157 0.063 0.002
[3,] 0.000 0.013 0.096
[4,] 0.003 0.073 0.069
[5,] 0.000 0.035 0.093
      [,1]  [,2]  [,3]
[1,] 0.010 0.031 0.010
[2,] 0.055 0.016 0.001
[3,] 0.000 0.003 0.038
[4,] 0.001 0.018 0.028
[5,] 0.000 0.009 0.037
[1] 0.051 0.071 0.042 0.047 0.046

　「関数による確率の計算」では、変数 $x$ のスカラとパラメータ $\boldsymbol{\lambda}$ のベクトルを指定して、複数のパラメータの計算を dpois() で行いました。
　今回は、変数 $\mathbf{x}$ を複製したマトリクスとパラメータ $\boldsymbol{\lambda}$ を複製したマトリクスを指定して、各変数 $x_n$ と各パラメータ $\lambda_k$ の組み合わせごとの計算を行います。dpois() に同じ形状のマトリクスやベクトルを指定することで、同じ位置の要素ごとに確率を計算できます。

　中間オブジェクトを作らずに、同様の計算を行います。

# クラスタごとのポアソン分布の確率を計算
cluster_prob_nk <- outer(X = x_n, Y = lambda_k, FUN = dpois)

#周辺確率を計算
prob_n <- (cluster_prob_nk %*% pi_k) |> 
  as.vector()
round(cluster_prob_nk, 3); round(prob_n, 3)

      [,1]  [,2]  [,3]
[1,] 0.029 0.125 0.025
[2,] 0.157 0.063 0.002
[3,] 0.000 0.013 0.096
[4,] 0.003 0.073 0.069
[5,] 0.000 0.035 0.093
[1] 0.051 0.071 0.042 0.047 0.046

　2つのベクトルの要素の全ての組み合わせに対する計算(処理)は、outer() で行えます。X, Y 引数にベクトル、FUN 引数に実行する関数を指定します。FUN 引数の関数の第1引数に X 引数のベクトル、第2引数に Y 引数のベクトルが複製されて渡されます。その他の引数を使う場合は、指定した関数の引数名で outer() に指定します。
　重み付け和の計算を行列演算で行います。

　作図など用に、データフレーム内で確率を計算します。

# クラスタごとの重み付け確率を計算
cluster_prob_df <- tidyr::expand_grid(
  n = 1:N, # データ番号
  k = 1:K, # クラスタ番号
) |> # 変数ごとにクラスタを複製
  dplyr::mutate(
    x      = x_n[n],      # 確率変数
    lambda = lambda_k[k], # パラメータ
    pi     = pi_k[k],     # 混合比率
    cluster_prob = dpois(x = x, lambda = lambda), # 各クラスタの確率
    weighted_prob = pi * cluster_prob             # 各クラスタの重み付け確率
  )
cluster_prob_df

# A tibble: 15 × 7
       n     k     x lambda    pi cluster_prob weighted_prob
   <int> <int> <dbl>  <dbl> <dbl>        <dbl>         <dbl>
 1     1     1    10    5.5  0.35    0.0285        0.00998  
 2     1     2    10   10    0.25    0.125         0.0313   
 3     1     3    10   16.8  0.4     0.0250        0.00998  
 4     2     1     6    5.5  0.35    0.157         0.0550   
 5     2     2     6   10    0.25    0.0631        0.0158   
 6     2     3     6   16.8  0.4     0.00158       0.000632 
 7     3     1    17    5.5  0.35    0.0000443     0.0000155
 8     3     2    17   10    0.25    0.0128        0.00319  
 9     3     3    17   16.8  0.4     0.0962        0.0385   
10     4     1    13    5.5  0.35    0.00277       0.000968 
11     4     2    13   10    0.25    0.0729        0.0182   
12     4     3    13   16.8  0.4     0.0690        0.0276   
13     5     1    15    5.5  0.35    0.000398      0.000139 
14     5     2    15   10    0.25    0.0347        0.00868  
15     5     3    15   16.8  0.4     0.0927        0.0371

　 $N$ 個の変数 $\mathbf{x}$ と $K$ 個のパラメータ $\boldsymbol{\lambda}$ の全ての組み合わせを expand_grid() で作成します。データ番号 $n = 1, \cdots, N$ とクラスタ番号 $k = 1, \dots, K$ の全ての組み合わせを作成して、 $\mathbf{x}, \boldsymbol{\lambda}, \boldsymbol{\pi}$ それぞれのベクトルから対応する値(要素)を取り出してデータフレームに格納します。各行が各組み合わせに対応します。
　mutate() で組み合わせ(行)ごとに重み付け確率を計算します。

　周辺確率を計算します。

# 周辺確率を計算
prob_df <- cluster_prob_df |> 
  dplyr::summarise(
    marginal_prob = sum(weighted_prob), # 周辺確率
    .by = c(n, x)
  )
prob_df

# A tibble: 15 × 7
       n     k     x lambda    pi cluster_prob weighted_prob
   <int> <int> <dbl>  <dbl> <dbl>        <dbl>         <dbl>
 1     1     1    10    5.5  0.35    0.0285        0.00998  
 2     1     2    10   10    0.25    0.125         0.0313   
 3     1     3    10   16.8  0.4     0.0250        0.00998  
 4     2     1     6    5.5  0.35    0.157         0.0550   
 5     2     2     6   10    0.25    0.0631        0.0158   
 6     2     3     6   16.8  0.4     0.00158       0.000632 
 7     3     1    17    5.5  0.35    0.0000443     0.0000155
 8     3     2    17   10    0.25    0.0128        0.00319  
 9     3     3    17   16.8  0.4     0.0962        0.0385   
10     4     1    13    5.5  0.35    0.00277       0.000968 
11     4     2    13   10    0.25    0.0729        0.0182   
12     4     3    13   16.8  0.4     0.0690        0.0276   
13     5     1    15    5.5  0.35    0.000398      0.000139 
14     5     2    15   10    0.25    0.0347        0.00868  
15     5     3    15   16.8  0.4     0.0927        0.0371

　summarise() のグループ化の引数 .by に変数 $x_n$ の列を指定して、変数ごとに重み付け確率の和を計算します。

　以上で、混合ポアソン分布の確率の計算を確認しました。

　この記事では、ポアソン分布の確率などの計算を確認しました。次の記事では、確率分布のグラフを作成します。

参考文献

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

作者:須山敦志
講談社

Amazon

おわりに

　確率分布シリーズの3回目の全編加筆修正中に新規で混合ポアソン分布編を追加しています。1記事に全部載せした方が分かりやすいだろう期と、目的ごとに記事を細分化した方が分かりやすいだろう期があって、今回の改修では1記事1目的の方針です。
　最初に書いたときの気持ち(目的や意図、方針など)がもう思い出せなくなりつつあり、確率分布の複雑な式を組んで理解しようだったのか、スクラッチ実装してパラメータの振る舞いを確認してから関数を使えるようにする流れにしたかったのか、確率の他に統計量やらの計算もまとめて1つの記事にするつもりだったのか、作図の解説の中で確率や確率密度の計算の説明も含めると流れが悪いので別の記事にしたのか、自分が分からなかったことをメモ的に並べただけなのか。
　今回の「確率分布の計算」の記事だと、「確率分布の作図」の記事で作図処理の解説に重きをおくために計算処理の解説をこっちで済ませておきたい気持ちです。ただそれだけだと書くことが少なく、1記事分の文量にするために水増し的に計算処理方法を追加してしまった気もします。
　基準は分かりませんがさすがに自分も中級者を名乗っていいだろうと思うようになったのですが、成長とともに初学者の気持ちが段々と想定できなくなってきて、何を書いてどこを解説すればいいのか悩むことが増えました。皆さんはいったい何を求めてこのブログに辿り着いたのでしょうか。

　さて投稿日の前日である2025年7月7日は、つばきファクトリーの3期(メジャーデビュー基準だと2期?)メンバーの加入4周年の日でした！

　また、Juice=Juiceの4期(あるいは6期?)メンバーの加入4周年の日でもありました！

　先日、新メンバーも加入して重ねてめでたい。

　えーっと、このブログでは何らかの記念日をよく取り上げるのですが、私はそういった日をまめに把握している性質ではなく、当日にSNSに流れてくるおめでとうの投稿を見て知り、それから書きかけの記事を完成させて文末にコメントを付け足して投稿する、という流れで運営されています。というわけで、間に合わなかったとか忘れていたとかいうよりもそもそも覚えておらず、そして今回はSNSを見る頻度が少ない日だったため知ったときには記念日を過ぎてしまっていました。
　開設からずっと記念日駆動執筆で運営されているブログなのですが、記念日に気付かないと更新されないというのは盲点でした。どうしたものか。

　とにもかくにも、おめでとうございます。私の日々の癒しの世代です。これからもよろしくお願いします。

【次の内容】

からっぽのしょこ

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

【R】混合ポアソン分布の計算