混合モデルによる潜在変数の特定(PyMC3)
1.背景
観測データの背後にある構造を把握することで有用な示唆を得られることがある。今回は簡単のために以下のような時系列データを考えてみる。下記はある人の心拍数の時間変位を想定して生成したデータである。今回、心拍数を決定づけるモードに平時、軽運動時、重運動時を想定し、観測されるデータから測定者が特定のタイミングでどのようなモードにあったのかを特定したい。混合モデルを使用することでモードの特定が可能になる。
2.混合モデル（Mixture Model）
混合モデル（Mixture Model）は、統計モデリングの手法の一つで、複数の確率分布の組み合わせでデータをモデル化する方法である。混合モデルは、複数の異なるデータ生成プロセスが同じデータセット内で共存している場合に適している。各データポイントは、異なる確率分布のいずれかから生成されたものと仮定される。
具体的な例を挙げると、ガウス混合モデル（Gaussian Mixture Model, GMM）がよく知られる。GMMは、複数のガウス分布（正規分布）を組み合わせてデータをモデル化する。各ガウス分布は、異なる部分集合のデータを表現し、それぞれが異なる平均値と分散を持つ。データポイントは、これらのガウス分布のどれかからランダムに選ばれると仮定される。
混合モデルは一般的に以下の要素により構成される。
成分（Component）: 各確率分布を成分と呼び、混合モデルではこれらの成分が複数存在する。各成分は異なる特性やパラメータ（平均、分散など）を持つ。
混合比率（Mixture Weights）: 各成分がデータセットにどの程度の割合で存在するかを表すパラメータ。これらの混合比率は、各成分がデータに対して寄与する程度を調整する。
確率密度関数（Probability Density Function, PDF）: 各成分に関連する確率密度関数が定義される。これは、ある成分が与えられたデータポイントを生成する確率を示す。
混合モデルは、クラスタリング、異常検知、特徴の抽出などのさまざまな応用分野で利用される。GMM以外にも、ポアソン混合モデルや指数混合モデルなど、様々な種類の混合モデルが存在する。データの複雑な構造をモデル化する際に有用な手法である。
3.PyMC3によるポアソン混合モデルの実装
人の心拍数の時間変位疑似データを生成し、PyMC3を使用した混合ポアソンモデルにより、特定のタイミングでどのような運動モードであったか特定してみる。
3.1.データの生成
下記のように疑似データを生成する。以下のコードでは、平時をλ=60、軽運動時をλ=90、重運動時をλ=110のポアソン分布にそれぞれ従うとし、状態遷移行列を定義し、各状態から別の状態へ状態遷移行列に基づいて遷移させ、その状態でのポアソン分布からのサンプリングによりデータを生成する。
import numpy as np import pymc3 as pm import theano.tensor as tt import matplotlib.pyplot as plt import seaborn as sns sns.set() # データ生成 np.random.seed(42) n_states = 3 n_observation = 200 # 遷移確率行列 transition_probs = np.array([[0.9, 0.05, 0.05], [0.05,0.9, 0.05], [0.05, 0.05, 0.9]]) # 隠れ状態の初期確率 initial_probs = np.array([1.0, 0.0, 0.0]) # 隠れ状態の生成 hidden_states = [np.random.choice(range(n_states), p=initial_probs)] for _ in range(n_observation - 1): hidden_states.append(np.random.choice(range(n_states), p=transition_probs[hidden_states[-1]])) # 観測モデルのパラメータ obs_means = np.array([60, 90 ,110]) # ポアソン分布の平均パラメータ # 観測データの生成 observed_data = np.random.poisson(obs_means[hidden_states])
以下が生成したデータと各時刻での運動モード（0:平時、1:弱運動、2:重運動）である。ベイズモデリングによる混合ポアソンモデルを用いることでこの特定の時刻の運動モードを推定する。
3.2.混合ポアソンモデルの実装
ベイズモデリングにより混合ポアソンモデルを実装する。特定の時刻の心拍数は３つの運動モードのうち、どれか１つの運動モードに属するとし、各モードはパラメータλに基づくポアソン分布に従うとしてモデル化する。また、それぞれのモードのポアソン分布のパラメータλの事前分布には一様分布、カテゴリカル分布の事前分布にはディレクレ分布を設定する。
K = 3 with pm.Model() as model: #混合比率 pi = pm.Dirichlet("pi", a=10*np.ones(K), shape=K) #クラスタ毎のポアソン分布のパラメータ mu = pm.Uniform("mu", lower=0, upper=60, shape=K) z = pm.Categorical("z", p=pi, shape=n_observation) x = pm.Poisson('obs', mu=mu[z], observed=observed_data)
モデル化ができたので、MCMC法を用いて事後分布のサンプリングをする。
with model: trace = pm.sample(2000, tune=1000, chains=4, random_seed=42)
サンプリング結果を可視化する。
_ = pm.traceplot(trace) _ = pm.plot_posterior(trace)
traceplotではサンプリングの分布が確認できる。
plot_posteriorでは特定の時刻で各モードに属する確率が可視化できる。
サマリーを確認する。
pm.summary(trace)
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
z[0] 0.000 0.016 0.000 0.000 0.000 0.000 8016.0 8016.0 1.0
z[1] 1.946 0.226 2.000 2.000 0.005 0.003 2378.0 2378.0 1.0
z[2] 2.000 0.016 2.000 2.000 0.000 0.000 8018.0 8000.0 1.0
z[3] 1.804 0.397 1.000 2.000 0.007 0.005 3655.0 3655.0 1.0
z[4] 2.000 0.011 2.000 2.000 0.000 0.000 8016.0 8000.0 1.0
... ... ... ... ... ... ... ... ... ...
pi[1] 0.369 0.034 0.302 0.431 0.000 0.000 4880.0 6001.0 1.0
pi[2] 0.276 0.032 0.218 0.336 0.000 0.000 6841.0 6050.0 1.0
mu[0] 58.796 0.959 56.966 60.578 0.011 0.008 7621.0 5323.0 1.0
mu[1] 98.230 1.404 95.687 100.914 0.023 0.016 3870.0 5212.0 1.0
mu[2] 137.769 1.887 134.301 141.397 0.031 0.022 3821.0 5276.0 1.0
最後に今回モデル化した混合ポアソンモデルによって推定される運動モードと実際の運動モードを比較してみる。 traceに観測データを持ちいて更新された事後データが格納されている。trace['z']に特定の時刻での各モードに属する確率が格納されているのでその中で最大となるindex（モード）を推定値とする。以下のコードにより特定の時刻でのモードの推定値を抽出できる。
mode_ = [] for i in range(0, n_observation): unique, freq = np.unique(trace['z'][:, i], return_counts=True) mode = unique[np.argmax(freq)] mode_.append(mode)
以下が実際の運動モード（左図）と推定した運動モード（右図）の比較である。隣接するモード間ではところどころ推定ミスはあるもののモードをよく推定できている。
