混合分布を考える

ポアソン混合モデルでは、ポアソン混合分布、つまり、いくつかのポアソン分布を合わせた分布を扱う。

下記に、3つのポアソン分布を混ぜたポアソン混合分布を描いてみた。

import numpy as np
from scipy.special import factorial
import matplotlib.pyplot as plt

x = np.linspace(start=0, stop=50, num=1000)
# ポアソン分布の尤度関数だけ作っておいた
uni_poisson = lambda x, lambda_: (lambda_ ** x) / factorial(x) * np.exp(-lambda_)
poisson1 = uni_poisson(x, 3)
poisson2 = uni_poisson(x, 10)
poisson3 = uni_poisson(x, 40)

# 合計が1になるような3つの重み(2次元標準シンプレックス)を作って足し合わせる
mixed_poisson = poisson1 * 0.2 + poisson2 * 0.3 + poisson3 * 0.5

plt.plot(x,mixed_poisson)

上図では、$ Poi(x|λ=3) $ と $ Poi(x|λ=10) $ と $ Poi(x|λ=40) $ が $ 0.2 : 0.3 : 0.5 $ の割合で混合している。

上図は、「3つのポアソン分布の尤度関数を重み付けして足し合わせた新しい尤度関数からのサンプルの分布」とみなすこともできる一方で、

「1番目のポアソン分布からのサンプルの分布と2番目のポアソン分布からのサンプルの分布と3番目のポアソン分布からのサンプルの分布を足し合わせた分布」ともみることができる。

後者を用いて考えると、順序だてて立式しやすい。

データ $ x $ がポアソン混合モデルから発生する概要を下記に記す。

１．クラスタ $ k $ のポアソン分布のパラメータ $ λ $ が、$ p(λ) $ からサンプルされる。

２．クラスタの混合割合から、x の発生するポアソン分布が選ばれる（クラスタ番号 $ k $ が確定する）。

３．クラスタ $ k $ のポアソン分布に従って $ x $ がサンプルされる。

１→２→３の順に説明を進めていこう。

１．クラスタ k のポアソン分布のパラメータ λ が、p(λ) から選ばれる。

クラスタの数だけ、ポアソン分布がある。ポアソン分布の数だけ、ポアソン分布のパラメータがあるわけで。

クラスタ $k$ のポアソン分布のパラメータを $λ_k$ と表すことにする。

$λ_k$ の事前分布としては、ポアソン分布の共役事前分布であるガンマ分布を仮定する。

各ポアソン分布(クラスタ)それぞれについて、所与データから学習したいので、クラスタの数 $ K $ だけ、ガンマ分布
$$
Gam(λ_k|a_k, b_k) = C_G(a, b) λ^{a - 1} e^{- b λ}
$$
$$ C_G(a, b) = \frac{b^a}{Γ(a)} $$
を用意する。

$ k = 1, 2, .., K $ についてガンマ分布 $ Gam(λ_k|a_k, b_k) $ から $ λ_k $ の実現値を得れば手順１は終了だ。

２．クラスタの混合割合から、x の発生するポアソン分布が選ばれる（クラスタ番号 k が確定する）。

クラスタの混合割合はデータを使って更新していきたいので、確率変数として扱おう。

k番目のクラスタの混合割合を $ \pi_k $ で表すことにする。

$$
\textbf{π} = (\pi_1, \pi_2, ,,, , \pi_K)
$$

ただし、クラスタの総数は $ K $ とする。

冒頭の例でいえば、

\begin{eqnarray*}
\textbf{π} 
&=& (\pi_1, \pi_2, \pi_3) \\
&=& (0.2, 0.3, 0.5) \\
\end{eqnarray*}

である。

ここで、 $ \bf{π} $ の事前分布を考えることにする。

要素を足し合わせて $ 1 $ になる $ K $ 次元ベクトルを生成したいなら、ディリクレ分布がおすすめ。(実装は最後に掲載)

$$
Dir(\textbf{π}|\textbf{α}) = C_D(\textbf{α}) \prod_{k=1}^{K} {\pi_k}^{\alpha_k - 1}
$$

ここで、$ \bf{α} $ は $ K $ 次元ベクトルのパラメータであり、各要素 $ \alpha_k $ は正の実数値である。

定数項は、

$$
C_D(\textbf{α}) = \frac{Γ(\sum_{k=1}^{K} \alpha_k)}{\prod_{k=1}^{K} Γ(\alpha_k)}
$$

であり、$ \pi_k $ ($ \textbf{π} $ のk番目の要素)の期待値は、

$$
\frac{\alpha_k}{\sum_{i=1}^{K} \alpha_i}
$$

である。

いまいちだな～という人のために、多項分布との類似性を少しばかり紹介したい。

多項分布：

\begin{eqnarray*}
Mult(\textbf{m}|\textbf{π}, M)
&=& M! \prod_{k=1}^{K} \frac{{\pi_k}^{m_k}}{m_k!} \\
&=& \frac{M!}{\prod_{k=1}^{K} m_k !} {\pi_k}^{m_k}
\end{eqnarray*}

である。多項分布は、$ K $ 個あるカテゴリの中から復元抽出を $ M $ 回行い、$ k $ 番目のカテゴリが $ m_k $ 回選択される確率を算出してくれる、高校数学でもおなじみの離散確率分布。

整数 $ n $ に対して、ガンマ関数は $ Γ(n) = (n - 1)! $ となることなどから、関数の形の類似性を感じてほしい。

多項分布は、例えば $ N=100 $ のヒストグラムの特定の一つのパターンが生起する確率を計算する。

この例を使えば、ディリクレ分布は、ヒストグラムの値(高さ)を正規化したうえで、特定パターンが生起する確率密度を計算する。連続確率分布なのでパターンは無数にあるのだけど...

多項分布が離散確率分布であり、ディリクレ分布が連続確率分布である点には留意しておきたい。

話が逸れすぎた。

てきとうな $ \textbf{α} $ で初期化されたディリクレ分布 $ Dir(\textbf{π}|\textbf{α}) $ から、クラスタ混合割合 $ \textbf{π} $ の実現値が得られる。

例えば3つのポアソン分布があるとして、$ \textbf{π} = (0.2, 0.3, 0.5) $ だったとする。

ここで、クラスタ1/2/3から選ぶ方法は？

カテゴリ分布を使えばいい。

カテゴリ分布：

$$
Cat(\textbf{s}|\textbf{π}) = \prod_{k=1}^{K} {\pi_k}^{s_k}
$$

カテゴリ分布は、$ K $ 個あるカテゴリから $ k $ 番目のカテゴリが選択される確率を算出する離散確率分布である。(実装は最後に掲載)

$ \textbf{s} $ は、$ K $ 次元ベクトルであり、一つの要素のみ$1$であり、そのほかの$K-1$個の要素は0である。

例えば、$K = 3$ であれば、$ \textbf{s} $ の実現値は $ (0,0,1), (0,1,0), (0,0,1) $ のいずれかである。

$ \textbf{π} $ は、$ K $ 次元ベクトルであり、$ k $ 番目の要素は、$ k $ 番目のクラスタが選択される確率を示している。

つまり、$ \textbf{π} $ の $ k $ 次元目の要素の値を $ \pi_k $ と表せば、 $ \sum_{k=1}^{K} \pi_k = 1 $ である。

カテゴリ分布は、ベルヌーイ分布

$$
Bern(x|\mu) = \mu^x (1 - \mu)^{1 - x}
$$

を多次元に拡張したものである。

後だしだが、カテゴリ分布に対して、ディリクレ分布は共役事前分布である。

つまり、カテゴリ分布のパラメータ $ \textbf{π} $ にディリクレ分布を仮定したとき、所与データから計算される $ \textbf{π} $ の事後分布はカテゴリ分布になる。

話がとんでもなく長くなってしまったが、、、

手順２．クラスタの混合割合から、x の発生するポアソン分布が選ばれる（クラスタ番号 k が確定する）。

では、

２－１．てきとうな $ \textbf{α} $ で初期化されたディリクレ分布 $ Dir(\textbf{π}|\textbf{α}) $ から、クラスタ混合割合 $ \textbf{π} $ の実現値を得る。

２－２．クラスタ混合割合 $ \textbf{π} $ をパラメータにもつカテゴリ分布 $ Cat(\textbf{s}|\textbf{π}) $ からカテゴリが選ばれる。

という流れで解析的にお話が進んでいくことになる、ということ。

３．クラスタ k のポアソン分布に従って x がサンプルされる。

手順１で各ポアソン分布の形状が決まり、手順２でデータの所属するクラスタが確定した。ここから $ x $ の実現値を得ればよい。

いったんここまでのまとめ

データが生成する過程をまとめよう。

１．選択されたクラスタ $k$ のパラメータ $ λ_k $ が $ Gam(λ_k|a_k, b_k) $ からサンプルされる。

２ー１．クラスタの混合割合 $ \textbf{π} $ が、$ Dir(\textbf{π}|\textbf{α}) = C_D(\textbf{α}) \prod_{k=1}^{K} {\pi_k}^{\alpha_k - 1} $からサンプルされる。

２－２．クラスタ番号を表すベクトル $ \textbf{s} $ が、$ Cat(\textbf{s}|\textbf{π}) = \prod_{k=1}^{K} {\pi_k}^{s_k} $ からサンプルされる。

３．$ Poi(x|λ_k) $ から $ x $ がサンプルされる。

ディリクレ分布、カテゴリ分布の実装

ディリクレ分布

$$
Dir(\textbf{π}|\textbf{α}) = C_D(\textbf{α}) \prod_{k=1}^{K} {\pi_k}^{\alpha_k - 1}
$$

ここで、$ \bf{α} $ は $ K $ 次元ベクトルのパラメータであり、各要素 $ \alpha_k $ は正の実数値である。

定数項は、

$$
C_D(\textbf{α}) = \frac{Γ(\sum_{k=1}^{K} \alpha_k)}{\prod_{k=1}^{K} Γ(\alpha_k)}
$$

であり、$ \pi_k $ ($ \textbf{π} $ のk番目の要素)の期待値は、

$$
\frac{\alpha_k}{\sum_{i=1}^{K} \alpha_i}
$$

である。

$ K = 3 $ (3次元)、$ \textbf{α} = (1, 2, 3), (5, 10, 15), (10, 20, 30) $ について描画してみる。

3次元のグラフを描画する方法はいくつかあるが、、、

今回は、$ a + b + c = 1 $ を満たす $ (a,b,c) $ を2次元の3角形で描画する　三角グラフ(ternary diagram)　を採用してみよう。

pythonのライブラリとして、python-ternaryを使用する。

# 三角グラフを描画するためのライブラリをインストールする
!pip install python-ternary

import numpy as np
from scipy.special import gamma

def Dirichlet(pi_vec, alpha_vec=[1, 2, 3]):
  # πベクトル、αベクトルをnumpyのndarrayとして扱う
  pi_vec = np.array(pi_vec)
  alpha_vec = np.array(alpha_vec)
  # Πは、各インデックスに対する値を格納したnp.ndarrayにnp.prod(・)を適用することで実装
  main = np.prod(pi_vec ** (alpha_vec - 1))
  C_D = gamma(np.sum(alpha_vec)) / np.prod(gamma(alpha_vec))
  return main * C_D

from functools import partial
# dirichlet_with_new_alpha = partial(Dirichlet, alpha_vec=[1,2,3])
import matplotlib.pyplot as plt
import matplotlib
# 三角グラフ描画用のライブラリをインポート
import ternary

alpha = [(1, 2, 3), (2, 4, 6), (4, 8, 12)]

# 描画用関数
def show_Dirichlet(alpha=(2, 2, 2)):
  # ただの分割数。
  scale = 100
  figure, tax = ternary.figure(scale=scale)
  figure.set_size_inches(10, 8)
  # sample_functionに引き渡されるタプルの中身は総和'1'の三つ組みタプル
  tax.heatmapf(func=partial(Dirichlet, alpha_vec=alpha), boundary=False, style="triangular", cmap='rainbow')
  tax.boundary(linewidth=2.0)
  tax.gridlines(color="black", multiple=10, linewidth=1)
#  tax.gridlines(color="blue", multiple=2, linewidth=0.5)

  # (r,p)の組と、期待値('Exp')を表示
  alpha = np.array(alpha)
  # np.array()に対しても、内法表記使えるみたいだね。便利。
  alpha_sum = alpha.sum()
  Exp = tuple([np.around((alpha_i / alpha_sum), decimals=2) for alpha_i in alpha])
  # 結局αはtupleに戻すんだけどね。表示用に。
  
  # Set Axis labels and Title
  fontsize = 25
  offset = 0.14
  tax.set_title(title=f'α = {tuple(alpha)}, Exp = {Exp}', loc='center', fontsize=fontsize)
  tax.bottom_axis_label("$\\alpha_1$", fontsize=fontsize, offset=offset)
  tax.right_axis_label("$\\alpha_2$", fontsize=fontsize, offset=offset)
  tax.left_axis_label("$ \\alpha_3 $", fontsize=fontsize, offset=offset)

  # multipy、デフォルトでは1みたい。
  tax.ticks(axis='lbr', linewidth=1, multiple=10)
  tax.clear_matplotlib_ticks()
  tax.get_axes().axis('off')
  tax.show()

# 描画するよ
show_Dirichlet(alpha[0])
show_Dirichlet(alpha[1])
show_Dirichlet(alpha[2])

期待値は

$$
\frac{\alpha_k}{\sum_{i=1}^{K} \alpha_i}
$$

となっている。およそ合っているように見える。

$ \alpha_3 $ の座標値が黄色い箇所の中心と必ずしも一致しないのは、分布の裾の引き方が等方的ではないためであろう。

カテゴリ分布

$$
Cat(\textbf{s}|\textbf{π}) = \prod_{k=1}^{K} {\pi_k}^{s_k}
$$

カテゴリ分布は、$ K $ 個あるカテゴリから $ k $ 番目のカテゴリが選択される確率を算出する離散確率分布である。

$ \textbf{s} $ は、$ K $ 次元ベクトルであり、一つの要素のみ$1$であり、そのほかの$K-1$個の要素は0である。

$ \textbf{π} $ は、$ K $ 次元ベクトルであり、$ k $ 番目の要素は、$ k $ 番目のクラスタが選択される確率を示している。

$ \textbf{π} $ の $ k $ 次元目の要素の値を $ \pi_k $ と表せば、 $ \sum_{k=1}^{K} \pi_k = 1 $ である。

\

結果は想像にたやすいと思うが、$ K = 3 $ (3次元)、$ \textbf{π} = (0.1, 0.3, 0.6), (0.3, 0.3, 0.4) $　として実装してみよう。

import numpy as np
from scipy.special import gamma

# s_vecは1 of K のベクトル
def Categorical(s_vec, pi_vec=[0.333, 0.333, 0.333]):
  s_vec = np.array(s_vec)
  pi_vec = np.array(pi_vec)
  # Πは、各インデックスに対する値を格納したnp.ndarrayにnp.prod(・)を適用することで実装
  # このCategorical関数では、s_vecをリストにした2次元リストが使われているので、axisを指定し、次元数の数と同じ次元の戻り値を返す。
  #main = np.prod(pi_vec ** (s_vec), axis=1)
  main = np.prod(pi_vec ** (s_vec), axis=0).flatten()
  return main

import matplotlib.pyplot as plt

pi = [(0.1, 0.3, 0.6), (0.3, 0.3, 0.4)]

# 描画用関数
def show_Categorical(pi):
  print(pi)
  # 地味にS_nのベクトルを生成するのって一番いい方法何？単位行列の分解はどうだ？
  # 3次元
  identity = np.eye(3)
  s_vec = np.split(identity, indices_or_sections=len(identity))
  s = np.arange(1, 4)
  y = Categorical(s_vec=s_vec, pi_vec=pi)
  # (r,p)の組と、期待値('Exp')を表示
  plt.title(label=f'π = {pi}', loc='center')
  plt.ylim([0, 1.0])
  #plt.plot(s,y)
  plt.bar(s,y)

# 描画するよ
plt.figure(figsize=(12,3))

plt.subplot(121)
show_Categorical(pi[0])

plt.subplot(122)
show_Categorical(pi[1])

お粗末さまでした。