Bayesian Analysis with PyMC 勉強ノート 7 混合モデル

Osvaldo Martin 著, 金子武久 訳
Pythonによるベイズ統計モデリング (Bayesian Analysis with Python)
2018-06-26
https://www.kyoritsu-pub.co.jp/book/b10003944.html

Christopher M. Bishop 著, 元田浩 栗田多喜夫 樋口知之 松本裕治 村田昇 監訳
パターン認識と機械学習 ベイズ理論による統計的予測 (Pattern Recognition and Machine Learning)
2012-01-20
https://www.maruzen-publishing.co.jp/item/?book_no=294524

Pythonによるベイズ統計モデリングを読んだので学習メモを整理します。学習中のノートなので正確ではないかもしれません。使われているコードをPyMC3からPyMC5に書き直しました。
また、せっかくなのでPRMLの本を参照しながら読むことにしました。式の文字の扱いなどはPRML風になるよう頭に収めようと思っています。

https://www.kyoritsu-pub.co.jp/book/b10003944.html

https://www.pymc.io/welcome.html

混合モデル

単純な分布では表現できない複雑な問題をモデリングすることができる。

• 有限混合モデル
• ゼロ過剰ポアソン分布・ゼロ過剰ポアソン回帰
• ロバストロジスティック回帰
• モデルベースクラスタリング
• 連続型混合モデル

次のライブラリを使用する。

import matplotlib.pyplot as plt
import math
import numpy as np
from scipy import stats
import pandas as pd
import seaborn as sns

import pymc as pm
import arviz as az

az.style.use("arviz-darkgrid")

1. ガウス混合モデル

過程や現象を正規分布など知れた正準分布(canonical distribution)で表わせない場合、分布の組み合わせを仮定するといい感じなことがある。
部分母集団がいくつかの組み合わせ(例えばデータ生成元の正規分布がいくつか存在する場合など)になっている場合はこれを直接的に表現できる。また、多峰性の分布を扱いやすく表現するにはガウス混合モデルGMM を使うといい。

頻度論的な使い方としてはカーネル密度推定KDE は馴染み深い。KDEではヒストグラムのかわりに、ヒストグラムが生成されていると仮定する連続的な分布を考えることで、ノンパラメトリックに推定する。具体的には各点に正規分布を割り当てて、全体で合計することで求められる。

混合モデルは、ある個数の部分母集団がある時に、どの分布かわからないけどそれに従うデータがある状況で、各データをfitする母集団のもとに割り当てるような時に使う。これは階層モデルで実現できる。

潜在変数(latent variable): 実際には観測できない確率変数。どの分布がどの観測データに対応するかを規定する。

clusters = 3
n_cluster = [90, 50, 75]
n_total = sum(n_cluster)
mu_real, sig_real = [9, 21, 35], [2, 2, 2]

y = np.random.normal(np.repeat(mu_real, n_cluster), np.repeat(sig_real, n_cluster))
fig, axes = plt.subplots(1, 2, figsize=(10,3))
sns.kdeplot(np.array(y), ax=axes[0])
sns.histplot(np.array(y), ax=axes[1])
plt.show()

このデータから要素の正規分布の数(k=3)を既知として、平均(mu)、分散(sig)とそれぞれのデータ点がどの正規分布に割り当てられるかを推定する。

with pm.Model() as model:
    p = pm.Dirichlet("p", a=np.ones(clusters))
    category = pm.Categorical("category", p=p, shape=n_total)
    mu = pm.Normal("mu", mu=[10, 20, 30], sigma=10)
    sig = pm.HalfCauchy("sig", 5)

    y_pred = pm.Normal("y_pred", mu=mu[category], sigma=sig, observed=y)
    trace = pm.sample(2000, chains=2, cores=1)
    
az.plot_trace(trace, var_names=["p", "mu", "sig"])

y_hat = pm.sample_posterior_predictive(trace, model=model)

fig, axes = plt.subplots(2, 1, figsize=(8,5))

az.plot_ppc(y_hat, num_pp_samples=20, ax=axes[0])
sns.kdeplot(y, c="b", ax=axes[1])
axes[1].set_xlim(0,45)

Dirichlet分布

β分布の一般化。3次元の場合、ベクトル $(p,q,r)$ を返すディリクレ分布は、制約 $r=1-(p+q)$ の上で返すベクトルの分布を定義できる。2つの背反な事象にはベルヌーイ分布で表現でき、その確率の事前分布にβ分布を用いた。これを一般化し、k個の背反な事象を表現するカテゴリカル分布に対し、その確率の事前分布にはディリクレ分布を使う。

[1] Bela A. Frigyik, Amol Kapila, and Maya R. Gupta
Introduction to the Dirichlet Distribution and Related Processes
UWEE Technical Report, UWEETR-2010-0006
http://mayagupta.org/publications/FrigyikKapilaGuptaIntroToDirichlet.pdf

[1]より、ディリクレ分布の例

PRML 2.2.1 より、関数は以下。カテゴリカル分布(多項分布)の共役事前分布であることが示されているが、確率密度関連に関してはどこかでまとめて復習しておきたい...

$$\text{Dir}(\bold{μ}|\bold{α}) = \frac{Γ(α_0)}{Γ(α_1) \dots Γ(α_K)} \prod_{k=0}^{K} μ_k^{α_k-1}$$

周辺化されたGMM

潜在変数、上記の例ではcategory のサンプリングは分布の無駄な探索を起こしてしまい遅くなる原因になる。潜在変数を迷惑変数として周辺化して消すことで、これを解消できる。

with pm.Model() as model:
    p = pm.Dirichlet("p", a=np.ones(clusters))
    mu = pm.Normal("mu", mu=[10, 20, 30], sigma=10, shape=clusters)
    sig = pm.HalfCauchy("sig", 5)

    y_pred = pm.NormalMixture("y_pred", w=p, mu=mu, sigma=sig, observed=y)
    trace = pm.sample(2000, chains=2, cores=1)
    
az.plot_trace(trace, var_names=["p", "mu", "sig"])

2. ポアソンモデルとカウントデータ

カウントデータ(count data): 離散的な非負数のデータを扱うにはポアソン分布を使う。

ポアソン分布

ある時間で発生する、独立事象の回数をの確率を記述するのにポアソン分布が使われる。確率質量関数(probability mass function; pmf)により、以下のように表せる。$λ$: 単位時間あたりの事象の生起期待値。

$$\text{Po}(k|λ) = \frac{λ^k}{k!} e^{-λ}$$

lambda_list = [0.5, 1.5, 3, 8]
k = np.arange(0, max(lambda_list)*3)
for l in lambda_list:
    y = stats.poisson(l).pmf(k)
    plt.plot(k, y, "-o", label=f"λ = {l:3.1f}")
plt.legend()

2項分布 $\text{Binomial}(k|n,p) = \frac{n!}{k!(n-k)!} p^k (1-p)^{n-k}$ を思い出すと、これは正規確率 $p$ の事象を $n$ 回独立に試行し、 $k$ 回生起した確率質量関数を示す。
ポアソン分布は、 $p = λ/n$ としたときの極限 $\lim_{n→∞} \text{Binomial}(k|n,λ/n)$ から導出される。

$$\begin{aligned} \text{Binomial}(k|n,λ/n) & = \frac{n!}{k!(n-k)!} (\frac{λ}{n})^k (1-\frac{λ}{n})^{n-k} \\ & = \frac{λ^k}{k!} (1-λ/n)^n (1-1/n)(1-2/n)(1-3/n) \dots (1-(k-1)/n)(1-λ/n) \\ & ↓ \lim_{n→∞} \\ \text{Po}(k|λ) & = \frac{λ^k}{k!} e^{-λ} \end{aligned}$$

ゼロ過剰ポアソンモデル

https://www.pymc.io/projects/docs/en/latest/api/distributions/generated/pymc.ZeroInflatedPoisson.html

ゼロ過剰ポアソンモデル(zero-inflated Poisson; ZIP)は、事後予測チェックなどを行う時に、データがポアソン分布に従っていると期待より過剰に0を観測してしまい、うまくデータにfitしないことがある問題を解決する。
ポアソン分布 $Ψ$ と過剰なゼロを生む確率 $(1-Ψ)$ により以下のように表せる。

$$\begin{aligned} \text{ZIP}(k|λ) & = Ψ \frac{λ^k}{k!} e^{-λ} \\ \text{ZIP}(0|λ) & = 1 - Ψ + Ψe^{-λ} \end{aligned}$$

この記事ではあまり深入りせず使ってみる。より専門的な、わかりやすい記事としては以下のNospareのものが良かった。

https://qiita.com/gen_nospare/items/20ecf330fbe2f26b7752

# ゼロが多いデータに対するZIPのfitting

n = 100
theta = 2.5
pi = 0.5
y = np.array([(np.random.random() > pi)*np.random.poisson(theta) for _ in range(n)])

with pm.Model() as model:
    psi = pm.Beta("psi", alpha=1, beta=1)
    lam = pm.Gamma("lambda", alpha=2, beta=0.1)
    y_pred = pm.ZeroInflatedPoisson("y_pred", psi=psi, mu=lam, observed=y)
    trace = pm.sample(2000, chains=2, cores=1)
    
az.plot_trace(trace, var_names=["psi", "lambda"])

fig, axes = plt.subplots(1, 2, figsize=(10,5))
y_sampled = pm.sample_posterior_predictive(trace, model=model)
az.plot_ppc(y_sampled, num_pp_samples=20, ax=axes[0])
sns.histplot(y, color="b", ax=axes[1])

ポアソン回帰・ゼロ過剰ポアソン回帰

線形モデルの一部としてポアソンやZIPを使うこともある。例えば、通常の実数範囲における一般化線形回帰モデルは活性化関数と線形モデルや誤差分布(たとえばガウスやT分布)を使って構築できた。分類のためのロジスティック回帰なども一般化線形モデルの枠組みで作れる。このアイディアと同じく、ポアソン分布を使って離散的な目標変数を回帰することができる。

# 訪問者の属性から、釣った魚の数を予測する

data = pd.read_csv("https://stats.oarc.ucla.edu/stat/data/fish.csv")
print(data)

x1, x2 = data["child"].values, data["camper"].values
y = data["count"].values

with pm.Model() as model:
    psi = pm.Beta("psi", alpha=1, beta=1)
    w = pm.Normal("w", mu=0, sigma=10, shape=3)
    lam = pm.math.exp(w[0] + x1*w[1] + x2*w[2])

    y_pred = pm.ZeroInflatedPoisson("y_pred", psi=psi, mu=lam, observed=y)
    trace = pm.sample(2000, chains=2, cores=1)

az.plot_trace(trace, var_names=["psi", "w"])
az.plot_autocorr(trace, var_names=["psi", "w"], combined=True)

fig, axes = plt.subplots(1, 3, figsize=(10,4))

for i, y_i in enumerate(y):
    if x2[i] == 0:
        axes[0].scatter(x1[i], y_i, color="c")
    else:
        axes[1].scatter(x1[i], y_i, color="m")  # キャンプ者
axes[0].set_ylim(0,40)
axes[1].set_ylim(0,40)

w_post_mean = trace.posterior["w"].mean(dim=["chain","draw"]).values
c = np.array([0,1,2,3])
axes[2].plot(np.exp(w_post_mean[0] + w_post_mean[1]*c + w_post_mean[2]*0), "-o", color="c")
axes[2].plot(np.exp(w_post_mean[0] + w_post_mean[1]*c + w_post_mean[2]*1), "-o", color="m")
plt.show()

output

     nofish  livebait  camper  persons  child        xb        zg  count
0         1         0       0        1      0 -0.896315  3.050405      0
1         0         1       1        1      0 -0.558345  1.746149      0
2         0         1       0        1      0 -0.401731  0.279939      0
3         0         1       1        2      1 -0.956298 -0.601526      0
4         0         1       0        1      0  0.436891  0.527709      1
..      ...       ...     ...      ...    ...       ...       ...    ...
248       1         1       1        3      2  1.374641 -2.595630      0
249       1         1       1        2      1  0.828834 -1.457115      0

ロバストロジスティック

目的変数に異常に大きな、あるいは小さな値がある場合でもそれに引っ張られずに分類するため、混合モデルを使える。
目的変数が $π$ でランダムに発生する(つまり $1-π$ で回帰モデルに依存する)という仮定をすることで、異常値かそうでないかを含めてモデルに組み込むことができる。

$$\begin{aligned} \text{Bernoulli}(x|p) & = p^x (1-p)^(1-x) \\ p & = 0.5π + (1-π) \text{σ}(\bold{w}^\top \bold{x}) \end{aligned}$$

$π = 1$ の時 $p = 0.5$ で一様分布になり、 $π = 2$ の時 $p$ はロジスティック回帰と等しくなる。

iris = sns.load_dataset("iris")
dataset = pd.concat([iris[iris["species"] == "versicolor"], iris[iris["species"] == "setosa"]])
y = np.concatenate([pd.Categorical(dataset["species"]).codes, np.ones(8)])
x = np.concatenate([dataset["sepal_length"].values, np.array([4.0,4.1,4.2,4.3,4.0,4.1,4.2,4.3])])

fig, axes = plt.subplots(1, 2, figsize=(10,4))

# 通常のロジスティック回帰
with pm.Model() as model:
    w = pm.Normal("w", mu=0, sigma=10, shape=2)
    mu = pm.Deterministic("mu", w[0] + pm.math.dot(x, w[1]))
    theta = pm.Deterministic("theta", 1/(1 + pm.math.exp(-mu)))
    bd = pm.Deterministic("bd", -w[0]/w[1])
    y_pred = pm.Bernoulli("y_pred", p=theta, observed=y)

    trace = pm.sample(2000, chains=2, cores=1)

idx = np.argsort(x)
post_theta = trace.posterior["theta"].values
post_theta_hdi = az.hdi(trace.posterior, hdi_prob=0.90)["theta"].values
post_bd = trace.posterior["bd"].values
post_bd_hdi = az.hdi(trace.posterior, hdi_prob=0.90)["bd"].values

axes[0].scatter(x, y, color="b")
axes[0].plot(x[idx], post_theta.mean(axis=(0,1))[idx], color="r")
axes[0].fill_between(x[idx], post_theta_hdi[idx][:,0], post_theta_hdi[idx][:,1], color="r", alpha=0.4)
axes[0].axvline(post_bd.mean(), ymax=1, color="r")
axes[0].fill_betweenx(y=[0,1], x1=post_bd_hdi[0], x2=post_bd_hdi[1], color="r", alpha=0.4)

# ロバストロジスティック
with pm.Model() as model:
    w = pm.Normal("w", mu=0, sigma=10, shape=2)
    mu = pm.Deterministic("mu", w[0] + pm.math.dot(x, w[1]))
    theta = pm.Deterministic("theta", 1/(1 + pm.math.exp(-mu)))
    bd = pm.Deterministic("bd", -w[0]/w[1])

    pi = pm.Beta("pi", alpha=1, beta=1)
    p = pi * 0.5 + (1 -pi) * theta
    y_pred = pm.Bernoulli("y_pred", p=p, observed=y)

    trace = pm.sample(2000, chains=2, cores=1)

idx = np.argsort(x)
post_theta = trace.posterior["theta"].values
post_theta_hdi = az.hdi(trace.posterior, hdi_prob=0.90)["theta"].values
post_bd = trace.posterior["bd"].values
post_bd_hdi = az.hdi(trace.posterior, hdi_prob=0.90)["bd"].values

axes[1].scatter(x, y, color="b")
axes[1].plot(x[idx], post_theta.mean(axis=(0,1))[idx], color="r")
axes[1].fill_between(x[idx], post_theta_hdi[idx][:,0], post_theta_hdi[idx][:,1], color="r", alpha=0.4)
axes[1].axvline(post_bd.mean(), ymax=1, color="r")
axes[1].fill_betweenx(y=[0,1], x1=post_bd_hdi[0], x2=post_bd_hdi[1], color="r", alpha=0.4)

print(trace.posterior["pi"].mean())

output

0.29522286  ← 混合確率 π

3. 連続型混合モデル

ロバストロジスティクス回帰は、ロジスティック回帰とランダム性の推測による連続型混合モデルの一種であることはわかりやすい。 $π$ の値がモデルの確率的な切り替えスイッチになっていて、混合係数(mixing coeffiicent)に相当する。
また、階層モデルも連続型混合モデルの一種で、各グループのパラメータが上位の連続分布から発生していると考えられる。

β2項分布

$$\text{BetaBinomial}(k|n,α,β) = \int_{0}^{1} \text{Binomial}(k|p,n) \text{Beta}(p|α,β) dp$$

これは $n$ 回の背反な試行(ベルヌーイ試行)において $k$ 回事象が生起する確率分布。
すなわち、 $k$ が観測される確率を表すために連続値 $p$ で期待値を取ることで、連続型混合モデルを形成している。
モデル内で期待値を取るパラメータは潜在変数のことを表している。

負の2項分布

連続型ガンマ-ポアソン混合(gamma-Poisson mixture)モデルの一つで、ガンマ分布による連続値を使ってポアソン確率 $p$ の期待値を取り、連続型混合モデルを形成する。
これにより、カウントデータによくある過剰分散 (データの分散がポアソン分布モデルの分散を超える)問題を解決できる。ポアソン分布の平均と分散はパラメータ1つを共有しているので、この部分をガンマ分布の情報を使って扱うことが良いらしい(?)

T分布

$$T_ν(y|μ,σ) = \int_{0}^{∞} \text{N}(y|μ,σ) \text{Invχ}^2(σ|ν) dν$$

$ν$ は $σ$ をサンプリングして発生する自由度(正規性パラメータ)。

感想

これまで具体的な混合モデルはGMMしか知らなかった(使っていたと思っていなかった)が、実際に式で見てみると階層モデルなどを始め知らないところで使っていた。というより現象のモデル(仮説)を考えるときに混合モデルになるのは自然なので、GMM以外を使うときも意識してモデリングに組み込めれば、説明性の観点からよりよい推論ができるようになるだろう。

最後の連続型混合モデルについては詳しく触れなかったのでまだ理解していない。PRMLだと下巻の14.5あたりに書かれている気がするので、読み追いついたらその時に詳しく見る予定。(下巻の難易度高い)