ガシンラーニング

マシンラーニング・統計

【Python実装】LDAのトピックをParticle Filter(SMC)で推論

 今回は、LDA(Latent Dirichlet Allocation)の逐次モンテカルロ法(Sequential Monte Calro)であるパーティクルフィルター(Particle Filter)によるトピック推論をPythonで実装しました。

コードは全てgithubに載せています。githubはこちら

Twitterフォローよろしくお願いいたします。twitterはこちら 

 

以下の書籍3.5章とこの書籍が参照している元論文を参考にしました。

Online Inference of Topics with Latent Dirichlet Allocation [Canini 2009]こちら

こちらの書籍はトピックモデルに限らずベイズモデリング推論の良書です。

トピックモデルによる統計的潜在意味解析 (自然言語処理シリーズ)

トピックモデルによる統計的潜在意味解析 (自然言語処理シリーズ)

初めに

膨大な量の文章が日々新たに流れてくる状況で、全ての文章の情報を保存して、学習の度に全データを読み込み、変分ベイズMCMCなどのバッチ学習をすることは非常に効率の悪くなる場面があります。そこで今回はミニバッチ学習でもなく、オンライン学習であるパーティクルフィルターによるメモリ効率の良いトピック推論を紹介したいと思います。

本記事を読むに当たって、LDAの崩壊型ギブズサンプリングによるトピック推論、つまり下の式の意味が分かることが望ましいです。導出過程は参考書籍に載ってます。

(文章\(d\)の\(i\)番目の単語が\(w_{d,i}\)=\(v\)である時の潜在トピック\(z_{d,i}\)の条件付き分布)$$p(z_{d,i}=k|w_{d,i}=v, {\bf{z}}^{\backslash{d,i}},{\bf{w}}^{\backslash{d,i}}, \alpha, \beta)$$

$$\propto \frac{{n_{k,v}}^{\backslash{d,i}}+\beta_v}{\sum_{{v}'}({n_{k,{{v}'}}}^{\backslash{d,i}}+\beta_{{v}'})}\frac{{n_{d,k}}^{\backslash{d,i}}+\alpha_k}{\sum_{{k}'}({n_{d,{k}'}}^{\backslash{d,i}}+\alpha_{{k}'})}$$

LDA(Latent Dirichlet Allocation)とは?

LDAの説明をしていませんでした。

LDAとは各文書に潜在的なトピック(意味)があると仮定した、文章の確率的生成モデルです。単語の相対的な出現頻度に注目しており、出現順序は無視しています。

f:id:gashin_learning:20191002115023p:plain

データの生成過程として、①「文章毎のカテゴリカル分布」からトピックが決まる。②「①で決まったトピックのカテゴリカル分布」から単語が決まる。という流れで、単語の数だけこの生成過程を辿ります。カテゴリカル分布のパラメータ\(\theta_d, \phi_k \)は各々ディリクレ分布から確率的に決まります。

Particle Filter パーティクルフィルターとは?

LDAに限らず、一般的な説明をします。

パーティクルフィルター(PariticleFilter)とは、逐次インポータンスサンプリング(Seqential Importance Sampling)リサンプリングを組み合わせた逐次ベイズフィルターのことで、潜在状態変数の事後分布\(p(z_{1:t}|x_{1:t})\)を多数の「パーティクル」と「重み」によって近似します。(\(z_{1:t}=\{z_1, z_2, z_3, \cdot \cdot \cdot ,z_t\}\)とする)

まずインポータンスサンプリングの説明から入ります。数式は参考書籍を基に記載し、適時補足します。

インポータンスサンプリング(Importance Sampling)

\(p(z_{1:t}|x_{1:t})\)が複雑な分布の場合、サンプルが容易な提案分布\(q(z_{1:t}|x_{1:t})\)からのサンプル\({{z}_{1:t}}^{(s)} \sim q(z_{1:t}|x_{1:t})\)(これがパーティクル)を用いて次のように近似します。

$$p(z_{1:t}|x_{1:t})=w(z_{1:t})q(z_{1:t}|x_{1:t}) \approx \tilde{ps}(z_{1:t}|x_{1:t})$$$$\tilde{ps}(z_{1:t}|x_{1:t})={\sum_{s=1}^{S} \frac{w({{z}_{1:t}}^{(s)})\delta({{z}_{1:t}}={{z}_{1:t}}^{(s)})}{\sum_{s=1}^{S} w({{z}_{1:t}}^{(s)})}}$$ここで$$w(z_{1:t})=\frac{p(z_{1:t}|x_{1:t})}{q(z_{1:t}|x_{1:t})}$$この\(w(z_{1:t})\)が重みと言われているものです。 

f:id:gashin_learning:20191002122057p:plain
 重みは正規化するため、\(p(z_{1:t}|x_{1:t})\)の数式が分からなければ、定数である\(p(x_{1:t})\)をかけて以下のようにしてもよいです。$$w(z_{1:t})=\frac{p(z_{1:t}|x_{1:t})}{q(z_{1:t}|x_{1:t})} \propto \frac{p(z_{1:t}, x_{1:t})}{q(z_{1:t}|x_{1:t})}$$同時分布にすればグラフィカルモデルに合わせて分解すれば求まります。

この\(\tilde{ps}(z_{1:t}|x_{1:t})\)を使うと、例えば\(p(z_{t+1}|x_{1:t})\)を求めたいときにも$$p(z_{t+1}|x_{1:t})=\int p(z_{t+1}|z_{1:t})p(z_{1:t}|x_{1:t})dz_{1:t}$$$$\approx \int p(z_{t+1}|z_{1:t})\tilde{ps}(z_{1:t}|x_{1:t})dz_{1:t}$$$$= {\sum_{s=1}^{S} \frac{w({{z}_{1:t}}^{(s)})p(z_{t+1}|{{z}_{1:t}}^{(s)})}{\sum_{s=1}^{S} w({{z}_{1:t}}^{(s)})}}$$ というように近似します。

 逐次インポータンスサンプリングでは、この\(\tilde{ps}(z_{1:t}|x_{1:t})\)から\({z_{1:t}}^{(s)}\)を一度にサンプルするのではなく、\({z_{c}}^{(s)} (c=1,2,3,,,t)\)を1つずつ逐次的にサンプルしていきます。

逐次インポータンスサンプリング(Seqential Importance Sampling)

 提案分布\(q(z_{1:t}|x_{1:t})\)は自由に選ぶことができるため、以下のような分解を仮定します。$$q(z_{1:t}|x_{1:t})=q(z_t|x_{1:t}, z_{1:{t-1}})q(z_{t-1}|x_{1:{t-1}}, z_{1:{t-2}})\cdot \cdot \cdot q(z_1|x_1, z_0)q(z_0)$$$$=q(z_0)\prod_{c}^{t}q(z_c|x_{1:c}, z_{1:{c-1}})$$つまり、$$q(z_{1:t}|x_{1:t}) = q(z_t|x_{1:t}, z_{1:{t-1}})q(z_{1:{t-1}}|x_{1:{t-1}})$$という提案分布の更新式を得ます。ここで数学的には\(q(z_{1:t}|x_{1:t}) = q(z_t|x_{1:t}, z_{1:{t-1}})q(z_{1:{t-1}}|x_{1:t})\)となりますが、\(q(z_{1:{t-1}}|x_{1:t}) = q(z_{1:{t-1}}|x_{1:{t-1}})\)と仮定します(提案分布なので自由に選べます。\(q(z_{1:t}|x_{1:t})\)にどのような構造を仮定するかは自由で、どのみち重みの更新式の計算で辻褄を合わせることになります)。提案分布の更新式は\(q(z_t|x_{1:t}, z_{1:{t-1}})\)になり、初期分布\(q(z_0)\)と\(t\)までの更新式の積として\(q(z_{1:t}|x_{1:t})\)を得ます。

重みも逐次的に更新します。重みの更新式は$$ \frac{ {w({z}_{1:t}}^{(s)}) } { {w({z}_{1:{t-1}}}^{(s)}) } = \frac{p({{z}_{1:t}}^{(s)}|x_{1:t})}{p({{z}_{1:{t-1}}}^{(s)}|x_{1:{t-1}})} \frac{q({{z}_{1:{t-1}}}^{(s)}|x_{1:{t-1}})}{q({{z}_{1:t}}^{(s)}|x_{1:t})}$$右辺第一項の分子分母にそれぞれ定数(\(x\)の確率分布)をかけて同時分布にします(どうせ後で正規化します)。右辺第二項は提案分布の更新式を使ってます。$$\propto \frac{p({{z}_{1:t}}^{(s)}, x_{1:t})}{p({{z}_{1:{t-1}}}^{(s)}, x_{1:{t-1}})} \frac{1}{q({z_t}^{(s)}|x_{1:t}, {z_{1:{t-1}}}^{(s)})}$$$$=\frac{p({z_t}^{(s)}, x_{t}|x_{1:{t-1}}, {z_{1:{t-1}}}^{(s)})}{q({z_t}^{(s)}|x_{1:t}, {z_{1:{t-1}}}^{(s)})} $$$$=\frac{p({x_t}|{z_t}^{(s)}, x_{1:{t-1}}, {z_{1:{t-1}}}^{(s)})p({z_t}^{(s)}|x_{1:{t-1}}, {z_{1:{t-1}}}^{(s)})}{q({z_t}^{(s)}|x_{1:t}, {z_{1:{t-1}}}^{(s)})} $$となります。

ちょっと脱線して、ブートストラップフィルター

ここで逐次的な提案分布を\(q({z_t}^{(s)}|x_{1:t}, {z_{1:{t-1}}}^{(s)})=p({z_t}^{(s)}|x_{1:{t-1}}, {z_{1:{t-1}}}^{(s)})\)とすれば、\({x_t}^{(s)}\)の情報が失われているので正確ではありませんが、重み更新式の一部分子分母が相殺されて\(p({x_t}|{z_t}^{(s)}, x_{1:{t-1}}, {z_{1:{t-1}}}^{(s)})\)だけが重み更新式として残ります。これが1次マルコフ状態空間モデルの場合、提案分布の更新式は「システムモデル」となり、重みの更新式は「観測モデル(尤度)」となりますね。これがブートストラップフィルターと呼ばれるパーティクルフィルターの特殊な場合になります。ブートストラップフィルターをパーティクルフィルターとして紹介している和書も多いです。(モンテカルロフィルターもこれと同じ?)

リサンプリング(Seqential Importance Sampling)

SISで逐次的に重みを更新していくと、サンプルサイズが増えるに従い、重みの分散が指数的に増加していきます。リサンプリングによって回避します。リサンプリングとは、重みの大きいパーティクルの生成と重みの小さいパーティクルの消滅を意味します。

f:id:gashin_learning:20191002124112p:plain

リサンプリングの後、元々の各パーティクルの重みは粒子の数として反映されたので、\(1/S\)にリセットされます。

リサンプリングではモンテカルロエラーを加えることになるため、毎回実施せずに\(EES\)と呼ばれる偏りの指標をもって、設定した閾値を下回った時に実施します。\(EES=\frac{1}{\sum_{s=1}^{S}(\bar{w}({z_{1:t}^{(s)})})^2}\)

リサンプリング手法も色々提案されていますが、今回のLDA推論では、元論文でも採用されているResidualResamplingにて行います。

(Liu and Chen 1998: "Sequential Monte Carlo Methods for Dynamic Systems"より)

f:id:gashin_learning:20191002124727p:plain
簡単に説明すると、「パーティクル数と正規化した重みとの積」を整数部分と少数部分に分けて、①整数部分の数だけ該当するパーティクルを複製する。②「元々のパーティクルの数」から「①で複製した後のパーティクルの合計数」を引いた数だけ、少数部分に比例した確率で抽出するといった手続きです。

 

一般的なパーティクルフィルターは以下の繰り返しです。

for t in range(data_num):

  1. 逐次的インポータンスサンプリングにより\({z_t}^{(s)}\)を\(S\)個サンプル
  2. 重み更新
  3. ESS計算。もしESSが閾値以下ならば、リサンプリング

LDAのトピックをParticleFilterで推論

パーティクルフィルターはマルコフ性のある状態空間モデルの潜在変数の推定に使われることが多いですが、今回LDAのパラメータ\(\theta_d, \phi_k \)を周辺化した場合に適用します。1次マルコフ性はないことに注意してください。(そのためパーティクルフィルタの説明でも1次マルコフ性を仮定せずに定式化していました。)

f:id:gashin_learning:20191002142933p:plain

LDAは時系列を仮定していないため文章の順番は適当に決めます。更新日時とかでいいんじゃないでしょうか。

前章までで、逐次的な提案分布の更新式は$$q(z_t|x_{1:t}, z_{1:{t-1}})$$重みの更新式は$$\frac{ {w({z}_{1:t}}^{(s)}) } { {w({z}_{1:{t-1}}}^{(s)}) } =\frac{p({x_t}|{z_t}^{(s)}, x_{1:{t-1}}, {z_{1:{t-1}}}^{(s)})p({z_t}^{(s)}|x_{1:{t-1}}, {z_{1:{t-1}}}^{(s)})}{q({z_t}^{(s)}|x_{1:t}, {z_{1:{t-1}}}^{(s)})} $$でした。

LDAのパーティクルフィルターでは提案分布を\(q(z_t|x_{1:t}, z_{1:{t-1}}) = p(z_t|x_{1:t}, z_{1:{t-1}})\)とし、これをLDAに合わせて書き直します。文章×単語の2次元行列を1次元に並べて$$p(z_{d,i}|w_{d,i}, {\bf{z}}^{(d, i-1)}, {\bf{w}}^{(d, i-1)}, \alpha, \beta)$$$$\propto p(w_{d,i}| z_{d,i}=k, {\bf{z}}^{(d, i-1)}, {\bf{w}}^{(d, i-1)}, \beta)p(z_{d,i}=k|\bf{z}^{(d, i-1)}, \alpha)$$$$\propto \frac{{n_{k,v}}^{(d,{i-1})}+\beta_v}{\sum_{v}({n_{k,v}}^{(d,{i-1})}+\beta_v)}\frac{{n_{d,k}}^{(d,{i-1})}+\alpha_k}{\sum_{k}({n_{d,k}}^{(d,{i-1})}+\alpha_k)}$$(この式は崩壊型ギブズサンプリングの式と同じ!ただし、\(x=w\)としています。\(w_{d,i}\)は文章dのi番目の単語。\(w^{(d, i)}=(w_{1:d}, w_{d,{1:i}})\)としています。)

重みの更新式は$$\frac{ w(z^{(d,i)(s)}) } { w(z^{(d,{i-1})(s)})} $$$$=\frac{ p(w_{d,i}| z_{d,i}, {\bf{z}}^{(d, i-1)}, {\bf{w}}^{(d, i-1)}, \beta)p(z_{d,i}|\bf{z}^{(d, i-1)}, \alpha)} {p(z_{d,i}|w_{d,i}, \bf{z}^{(d, i-1)}, \bf{w}^{(d, i-1)}, \alpha, \beta)} $$$$=\sum_{k=1}^{K}p(w_{d,i}| z_{d,i}=k, {\bf{z}}^{(d, i-1)}, {\bf{w}}^{(d, i-1)}, \beta)p(z_{d,i}=k|{\bf{z}}^{(d, i-1)}, \alpha)$$$$= \sum_{k=1}^{K} \frac{{n_{k,v}}^{(d,{i-1})}+\beta_v}{\sum_{v}({n_{k,v}}^{(d,{i-1})}+\beta_v)}\frac{{n_{d,k}}^{(d,{i-1})}+\alpha_k}{\sum_{k}({n_{d,k}}^{(d,{i-1})}+\alpha_k)}$$これも崩壊型ギブズサンプリングの式を\(z_{d,i}\)で周辺化しただけであるため計算容易です。

冒頭でパーティクルフィルターによる逐次学習によって、オンラインで更新された分のデータだけを使ってトピックを推論することができますと書きましたが、実は単語毎のトピック頻度\({n_{k,v}}^{(d,{i-1})}\)は保持する必要があります。ただこれは「ユニーク単語数×トピック」の行列で表現できるので文章が増えてもそんなに変わりにくいです。代わりに\({n_{d,k}}^{(d,{i-1})}\)は保持しなくてよいのが利点です。こちらは「文章数×トピック」なので文章数が容量に直結するからです。

リサンプリングについては、書籍ではリサンプリングを実施しないと書いていますが、元論文では実施しているためResidual Resamplingを実装しました。

元論文では、この後「若返り(rejuvenate)」というステップがあります。これは過去のトピックのサンプルを今時点までの全情報を使い、崩壊型ギブズサンプリングで再度サンプルし直すという手続きです。パーティクルの多様性を保つことが可能と説明していますが、後続論文でパフォーマンスとしてはあまり意味が無いとなっています。文章の潜在変数情報を記憶しておく必要があり、メモリ削減の利点が薄れるため実装には入れていません。 若返りたい人は崩壊型ギブズサンプリングの実装はあるため、文章のトピックを保持するように変えれば付け足し可能です。

実装

パーティクルフィルターはNumpyで実装しました。データ前処理&評価可視化にはscikit-learnを使っています。

パーティクルフィルタのアルゴリズムだけオンライン処理しています。文章の読み込み&ベクトル化の箇所はオンラインにしていないので、ここは用途に合わせてください。

コードは全てGithubに載せています。Githubはこちら

 SIS
### Sequential Importance Sampling

# p(z_k| z_{1:k-1})
n_d_alpha_cond_d_particles = np.sum(n_dk_alpha_cond_d_particles, axis=0)
p_z = n_dk_alpha_cond_d_particles /n_d_alpha_cond_d_particles

## p(w_k| z_{1:k}, w_{1:k-1})
n_vk_beta_cond_k_particles = n_vk_beta_particles[word_id, :, :]
n_k_beta_cond_k_particles = np.sum(n_vk_beta_cond_k_particles, axis=0)
p_w = n_vk_beta_cond_k_particles /n_k_beta_cond_k_particles

# p(z_k| z_{1:k-1})p(w_k| z_{1:k}, w_{1:k-1})
sample_p = p_z * p_w
sample_p_sum = np.sum(sample_p, axis=0)
sample_p /=sample_p_sum

# sampling from multinomial distribution
random_sampling = np.random.uniform(0,1, size=particle_num)
sample_p_cumsum = np.cumsum(sample_p, axis=0)
topic_new = multinomial_particles_numba(topic_new, random_sampling, sample_p_cumsum, particle_num)
topic_particles[idx] = topic_new

# update parameters
n_dk_alpha_cond_d_particles[topic_new, np.arange(100)] += 1.0
n_vk_beta_particles[word_id, topic_new, np.arange(100)] += 1.0

# update weight
w*= sample_p_sum
w /= np.sum(w)

# ESS
ESS = 1/np.sum(w**2)
resampling
### Residual Resampling 
K = particle_num * w
K_int = np.trunc(K)
redisual_p = K - K_int
redisual_p/=np.sum(redisual_p)

# main resample K_int = K.astype('int')
# residual resample
M = particle_num - np.sum(K_int)
residual_resampling_array = np.random.multinomial(n=M, pvals=redisual_p)
# resample process resampling_idx = K_int + residual_resampling_array particle_idx = np.repeat(np.arange(particle_num), resampling_idx) topic_particles = topic_particles[:,particle_idx]

# parameter update n_vk_beta_particles = n_vk_beta_particles[:,:,particle_idx] n_dk_alpha_cond_d_particles = n_dk_alpha_cond_d_particles[:,particle_idx]
# weight initizalization w = np.ones(particle_num)/particle_num

実験

データセット

データセットはsklearnのfetch_20newsgroupsを使用しています。簡単のため、このうち'rec.sport.baseball', 'talk.religion.misc','comp.graphics', 'sci.space'のみ使用しました。そのためトピック数も4に設定しています。

出現頻度min_df=0.005, max_df=0.1で、かつ名詞/固有名詞のみを抽出しています。

推論について

実際には、最初のデータの20%を使って崩壊型ギブズサンプリングでパラメーターを初期化します。このステップは精度の面で重要とされています。

ワードクラウド 

推論されたトピック毎の頻出ワードを列挙します。

f:id:gashin_learning:20191103162031p:plain

topic0は宇宙。topic1はパソコン。topic2は野球。topic3は宗教と元のニュースと同じトピックを抽出できています。

文章のラベルプロット

元データにどのニュース記事かのラベルがついているため(推論には使用していない)、t-SNEで2次元にして可視化して、ラベルとその文章に任意のワードが追加された時の予測トピックとの類似度をみました。

f:id:gashin_learning:20191103165515p:plain
色がラベル(どのニュース記事か)。座標が推論トピックの事後確率(4次元)を2次元に圧縮した時の位置。

所々ミスしていますが、一応色のかたまりが見てとれます。
 

高度化に向けて

元論文で提案されているパーティクル毎のトピックを木構造で格納することでメモリ効率を上げる手法は今回実装していないので残課題です。

 最後に

パーティクルフィルターの利点はデータに対して、一回しかサンプル(一回で粒子数サンプルする)しないため、推論に使用し終わった後はデータを保持しておく必要がないところです。

ただ、メモリに乗るくらいの文章量ならバッチ学習で全然構いません。その方が全てのデータを加味して推論できるため、一般的には精度も良いと思います。

場面に合わせた推論手法を選ぶと良いと思います。

Appendix ~ NumbaによるPython高速化~

今回、一部Numbaで高速化しました。

NumbaとはJITコンパイラで、既存のPythonコードに少し手を加えるだけで高速化できます。コードを事前にコンパイルすることで、ForループのNumpy処理が爆速になります。https://github.com/numba/numba

最も簡単な使い方は以下のようにnumbaをimportして、関数にデコレータをつけるだけです。

# using Numba 
@jit(nopython=True)
def collapsed_gibbs_sampling():
~

# Not using Numba
def collapsed_gibbs_sampling():
~

例えば崩壊型ギブズサンプリングをNumbaで高速化した場合をNumba使わない場合とで速度比較すると100倍以上速くなっています。

f:id:gashin_learning:20191103173909p:plain

ただ辞書やpandasの処理などは変換が上手くいかないようです。あくまでNumpyの高速化で、Numpyでもサポートされていない関数が結構あります。。。axis系だめ。。。。そのためパーティクルフィルタには使いませんでした。

[サポート関数一覧]

http://numba.pydata.org/numba-doc/latest/reference/numpysupported.html

今後、少しずつ使っていきたいです。

 

Twitterフォローよろしくお願いいたします。twitterはこちら 

【Python実装】ノンパラベイズ3次元無限関係モデル(3D-IRM)をギブスサンプリング(MCMC)で推論

今回は、書籍「続・わかりやすいパターン認識」の13章で紹介されている無限関係モデル(Infinite Relational Model)ギブズサンプリング(MCMCによる推論を、3次元にカスタマイズした3D-IRM(勝手に名前)をPythonで実装します。

モデルと推論方法に関しては、書籍「続・わかりやすいパターン認識」の13章を参考にしています。詳しくはこちらをご参照ください。

続・わかりやすいパターン認識―教師なし学習入門―

続・わかりやすいパターン認識―教師なし学習入門―

今回のコードを全てgithubに載せています。遊べるようにnotebookもつけてます。githubはこちら

Twitterフォローよろしくお願いいたします! twitterはこちら

 無限関係モデルとは?

詳しくは書籍「続・わかりやすいパターン認識」の13章を参考にしてください。ここでは簡単に説明します。

無限関係モデル(Infinite Relational Model:IRM)とは、異種オブジェクトを同時にクラスタリングする共クラスタリング手法の一つで、ノンパラメトリックベイズに基づきます。

ざっくり言うと、データの値は「1」か「0」で、軸毎にカチャカチャ入れ替えて、「1」が比較的多いグループと「0」が比較的多いグループにクラスタリングするイメージです。

下図が3D-IRMによるクラスタ結果の例です。「1」のみをプロットしています。ルービックキューブみたいです。

f:id:gashin_learning:20190906180832p:plain

ポイントは「X軸を動かすときはY, Z軸は固定」「Y軸を動かすときはX, Z軸は固定」と動かし方に制限があるところです。 

元のトイデータを4×4×4のキューブで全て統一していますが、データ次第ではクラスタ毎にキューブの大きさがバラバラになります。むしろ大半がそうです。(同軸上の長さは固定)

f:id:gashin_learning:20190906205704p:plain
ノンパラメトリックベイズなので、クラスタ数は予め決めず、データから自動で推論してくれます。書籍では、中華料理店過程(Chinese Restaurant Process、略してCRP)に従って独立に各軸がクラスタリングされると仮定しています。

仮定するデータの生成過程としては数式で表すと(各定義は書籍)

$${\bf s}^{1}| \alpha \sim C R P(\alpha)$$$${\bf s}^{2}| \alpha \sim C R P(\alpha) $$

$${\bf s}^{3}| \alpha \sim C R P(\alpha) $$$$ \theta({s_x}^{1}, {s_y}^{2}, {s_z}^{3})| a,b \sim Be(a,b)$$

$$R_{xyz}|{s_x}^{1} = {w_i}^{1}, {s_y}^{2} = {w_j}^{2}, {s_z}^{3} = {w_k}^{3} \Theta \sim Bern(R_{xyz} ; \theta_{ijk})$$

CRPによってクラスタ番号を表す\({\bf s}^{1}, {\bf s}^{2}, {\bf s}^{3}\)が割り当てられ、3軸の各々のクラスタの直積からなる各クラスタ\((i, j, k)\)毎にパラメータ\(\theta_{ijk}\)がベータ分布から別々に決まり、1セルずつ\(R_{xyz}\)(「1」か「0」)がベルヌーイ分布から生成されるといった具合です。

何に使えるのか?

例えば、(x, y, z)軸に(店舗、顧客、商品)をセットして、購入したを「1」、購入しないを「0」とすると、各クラスターから「この(地域A)では、(顧客層B)に対して、(商品カテゴリーC)が売れやすい、売れにくい」などが分かります。(x, y, z)に何をセットするかは自由です。応用幅は広いと思います。 

 

モデルに関しては以上で、推論方法について書籍では、崩壊型ギブスサンプリングで数式と手続きが書かれているのですが、今回は自分の勉強のために普通のギブスサンプリングに直して実装しました。

ギブズサンプリングとは?

MCMCの一手法で、各パラメーター毎にそれ以外のパラメーターで条件付けした分布から、順にサンプリングしていく手法です。

パラメーター毎に条件付き分布を手計算する必要があり、stanでよく使われる自動微分すればOKなHMC(Hamiltonian Monte Carlo)とは違い、自動化が難しい手法になっています。

ギブズサンプリングは、推移核を工夫したメトロポリス・ヘイスティング法(MH法)とも解釈ができ、MH法で言うところの採択率が1になり、サンプルが棄却されないことで知られています。1パラメーターずつサンプルしていくため、各軸平行に直角にサンプル点が動いていくイメージです。そのためパラメーター間の相関が強い分布に対しては、サンプリング数が多く必要になります。また、パラメーターの条件付き分布がいつも簡単にサンプリングできる分布になるとは限らないので要注意です。

f:id:gashin_learning:20190831162434p:plain

クラスタ番号\({\bf s}\)の条件付き分布が中華料理店過程にちょうど対応します。

ギブズサンプリングによってサンプリングした\({\bf s}^{1}, {\bf s}^{2}, \Theta\)の内、対数事後分布\(logP({\bf s}^{1}, {\bf s}^{2}, \Theta\|R)\)を最も大きくする\({\bf s}^{1}, {\bf s}^{2}, \Theta\)を最終結果として採用します。

実装

ほとんどNumpyとScipyで実装しました。データも自分で作成しています。出来るだけfor文を使わずにNumpyで書くようにしたので可読性は低いです(Numpyのせいにする)。それでも次元を増やしすぎると極端に遅くなります。

3D-IRMによる共クラスタリングアルゴリズムの箇所です。 

# gibbs sampling
def predict_S(R, alpha,a,b, iter_num=500, reset_iter_num=100):

    X, Y, Z = R.shape

    # set first values  ##########################
    sx = CRP(alpha=alpha, sample_num=X)
    sy = CRP(alpha=alpha, sample_num=Y)
    sz = CRP(alpha=alpha, sample_num=Z)
    theta = posterier_theta(sx, sy, sz, R, a, b)
    ##############################################

    max_v = -np.inf
    # to recycle 'def s_update'
    R_transpose_y = R.transpose((1,2,0))
    R_transpose_z = R.transpose((2,0,1))

    # gibbs sampling
    for t in range(iter_num):
        print("\r calculating... t={}".format(t), end="")
        sx, theta = s_update(sx, sy, sz, theta, R, a, b, alpha, axis=0)
        sy, theta = s_update(sy, sz, sx, theta, R_transpose_y, a, b, alpha, axis=1)
        sz, theta = s_update(sz, sx, sy, theta, R_transpose_z, a, b, alpha, axis=2)

        log_p_sx = np.log(Ewens_sampling_formula(sx, alpha))
        log_p_sy = np.log(Ewens_sampling_formula(sy, alpha))
        log_p_sz = np.log(Ewens_sampling_formula(sz, alpha))
        log_p_theta = np.sum(st.beta.logpdf(theta, a,b))

        log_p_R_theta = log_R_theta_probability(sx, sy, sz,  R, theta)
        #log_p_R_ijk = log_R_probability(sx, sy, sz,  R, a, b)

        # logP(sx, sy, sz, theta| R)
        v = log_p_sx + log_p_sy + log_p_sz + log_p_theta + log_p_R_theta

        # update if over max
        if v > max_v:
            max_v = v
            max_sx = sx
            max_sy = sy
            max_sz = sz
            max_theta = theta
            print("  update S and theta : logP(sx, sy, sz, theta| R) = ", v)

        # to prevent getting stuck local minima, reset S and theta
        if t%reset_iter_num==0:
            sx = CRP(alpha=alpha, sample_num=X)
            sy = CRP(alpha=alpha, sample_num=Y)
            sz = CRP(alpha=alpha, sample_num=Z)
            theta = posterier_theta(sx, sy, sz, R, a, b)

    return max_sx, max_sy, max_sz, max_theta

早い段階で更新されなくなり局所解にはまっている印象を受けたので数十サンプル毎にパラメーターを初期化して、広範囲探索するようにしました。

具体的なギブズサンプリングの部分です。

# update s
def s_update(s1, s2, s3, theta, R, a, b, alpha, axis):

    if axis==1:
        theta = theta.transpose((1,2,0))
    elif axis==2:
        theta = theta.transpose((2,0,1))
# sort orderby s2,s3 for easy calculation sorted_s2_index = s2.argsort() sorted_s3_index = s3.argsort() sorted_s2 = s2[sorted_s2_index] sorted_s3 = s3[sorted_s3_index] R_sorted = R[:,sorted_s2_index,:][:,:,sorted_s3_index] for idx in range(len(s1)): # remove s1_x for gibbs sampling s1_delete = s1[idx] s1_left = np.delete(s1, idx) theta_left = theta[np.unique(s1_left),:, :] # if category_num is decreased, fill empty category number s1_left = reset_s_number(s1_left) # count n_ij num_n_ijk_left = count_n_ijk(s1_left, s2, s3) # log_p(s1_k | s1_left) by Dirichlet Process n_i = np.add.reduce(num_n_ijk_left, axis=(1,2)) n_i = np.append(n_i, alpha) ln_p_s1_idx_s1_left = np.log(n_i/(np.add.reduce(num_n_ijk_left,axis=(0,1,2)) + alpha)) # log_p(R| s1_left, s2, s3) R_idx_sorted =R_sorted[idx,:,:] R_ijk_1, R_ijk_0 = count_one_zero_2D(R_idx_sorted, sorted_s2, sorted_s3) ln_p_R_xyz_new = np.sum(betaln(R_ijk_1+a, R_ijk_0 +b) - betaln(a,b)) ln_p_R_xyz_exist= np.sum(R_ijk_1 * np.log(theta_left), axis=(1,2))+ np.sum(R_ijk_0 * np.log(1-theta_left),axis=(1,2)) # Ratio for choosing new s1_x '+100' is for preventing underflow p_s1_idx = np.exp(ln_p_s1_idx_s1_left + np.append(ln_p_R_xyz_exist, ln_p_R_xyz_new)+100) p_s1_idx/=np.sum(p_s1_idx) s_new = np.argmax(np.random.multinomial(n=1, pvals=p_s1_idx)) # new s1 updated s1 = np.insert(s1_left, idx, s_new) # update theta theta = posterier_theta(s1, s2, s3, R, a,b) if axis==1: theta = theta.transpose((2,0,1)) elif axis==2: theta = theta.transpose((1,2,0)) return s1, theta

 今回、ギブズサンプリングを使用していますが、ベルヌーイ分布のパラメーターだけはサンプリングせず、条件付き分布の平均値をサンプルとみなして使っています。特に理由はないです。(この手法もなんか名前があるらしい)サンプリングしてもOKです。また軸毎に\(s\)をサンプルする毎にベルヌーイ分布のパラメーターもサンプル。\(s_x\)→\(\Theta\)→\(s_y\)→\(\Theta\)→\(s_z\)→\(\Theta\)と更新しています。

 

今回のコードを全てgithubに載せています。遊べるようにnotebookもつけてます。githubはこちら

動かしてみた

まずは、見やすさのためベルヌーイ分布のパラメーターを1か0にした超理想的なデータに適用しました。クラスタリング手法なのでデータによって、局所解が多く存在します。 f:id:gashin_learning:20190906180848p:plain

 次にノイジーなデータにも試しました。 パラメーターが0.8と0.2と0の3種類のベルヌーイ分布で生成しました。

f:id:gashin_learning:20190906212032p:plain

クラスタ数(4,4,4)で推定できました。X軸で切ってクラスタ毎のパラメータ\(\Theta\)を可視化しました。0.2や0.8付近ででまとまっています。

f:id:gashin_learning:20190906210847p:plain

今回のコードを全てgithubに載せています。遊べるようにnotebookもつけてます。githubはこちら

Twitterフォローよろしくお願いいたします! twitterはこちら

【Python実装】HMM(隠れマルコフモデル)を構造化変分ベイズで推論

今回は、HMM(隠れマルコフモデル構造化変分ベイズ推論Pythonで実装します。

モデルと推論方法に関しては、書籍「ベイズ推論による機械学習(須山)」の5-3章を参考にしています。詳しい途中計算はこちらをご参照ください。  

機械学習スタートアップシリーズ ベイズ推論による機械学習入門 (KS情報科学専門書)

機械学習スタートアップシリーズ ベイズ推論による機械学習入門 (KS情報科学専門書)

  • 作者: 須山敦志,杉山将
  • 出版社/メーカー: 講談社
  • 発売日: 2017/10/21
  • メディア: 単行本(ソフトカバー)

今回のコードをgithubに載せています。 githubはこちら

 

Twitterフォローよろしくお願いいたします! twitterはこちら 

HMM(隠れマルコフモデル)とは?

HMMとは、観測データの背後にある潜在変数\({\bf S}=\{s_1,s_2,...,s_N\}\)の離散的な系列変化を仮定したモデルです。時系列データに対するモデリングとして広く利用されています。

f:id:gashin_learning:20190825011936p:plain

観測モデル\({\bf X}\)と潜在変数の状態系列モデル\({\bf S}\)から成り立ちます。

今回は簡単のため、\(s_n\)が1時点前の\(s_{n-1}\)のみに依存する1次マルコフモデルを仮定します。

状態系列モデル

$$p(s_n|s_{n-1}, {\bf A})=\prod_{i = 1}^{K} Cat({\bf s}_n|{\bf A}_{:,i})^{s_{n-1}, i}$$$$p(s_1|\boldsymbol\pi)=Cat(s_1|\boldsymbol\pi)$$初期確率\(\boldsymbol\pi\)と遷移確率行列の列\({\bf A}_{:,i}\)にも各々ディリクレ分布で事前確率を与えます。

実務では\({\bf A}_{:,i}\)について、ドメイン知識を入れ込める場合が多いです。十分な学習に必要なデータ数の削減や精度の安定化ができることがあり、ベイズ学習の強みの一つといえます。

観測モデル

今回は書籍に則り、ポアソン観測モデルによる隠れマルコフモデルを仮定します。$$p(x_n|s_n)=\prod_{k = 1}^{K} Poi(x_n|\lambda_k)^{s_n, k}$$\(\lambda_k\)にも共役事前分布のガンマ分布で事前確率を与えます。観測モデルは、タスクに応じて自由に設定します。 

観測変数\({\bf x}\)のみを使用して、未知の潜在変数\({\bf S}\)と各パラメーター\(\boldsymbol\lambda, {\bf A}, \boldsymbol\pi\)の事後分布を構造化変分ベイズ推論します。

構造化変分ベイズとは? 

まず変分ベイズとは、複雑な事後分布を簡単な分布で近似するために2分布間のKLダイバージェンスを最小化する、またはELBO(変分下限)を最大化する推論方法です。最適化問題に帰着するため、MCMCと比べて計算が高速です。

変分ベイズでは計算を簡易化するため、推論したい潜在変数・パラメーターを全て独立と仮定する「平均場近似」による「完全分解変分ベイズ」が一般的に利用されています。

しかしそれでは、隣り合う潜在変数同士が相関を持つというHMMの仮定に反することになるため、今回は潜在変数同士は独立に分解せずに同時分布として扱う「構造化変分ベイズ」にて推論を実施します。

  • 完全分解変分ベイズ$$p({\bf S}, \boldsymbol\lambda, {\bf A}, \boldsymbol\pi)  \approx \{\prod_{n = 1}^{N} q(s_n)\} q(\boldsymbol\lambda, {\bf A}, \boldsymbol\pi)$$
  • 構造化変分ベイズ(今回実施。\(q({\bf S})\)が同時分布)$$p({\bf S}, \boldsymbol\lambda, {\bf A}, \boldsymbol\pi)  \approx q({\bf S})q(\boldsymbol\lambda, {\bf A}, \boldsymbol\pi)$$

フォーワード・バックワードアルゴリズム

\(q({\bf S})\)の同時分布を求めるといっても、\(q(\boldsymbol\lambda, {\bf A}, \boldsymbol\pi)\)の計算に\(s_n\)と \(s_{n-1} \rightarrow s_n\)の期待値が必要で、\(q(s_n), q(s_n, s_{n-1}),\)を使って算出するのですが、構造化変分ベイズの計算過程の中で\(q(s_n), q(s_{n-1},s_n)\)を個別に計算できるわけではないため、\(q({\bf S})\)を求めた後に周辺化して求める必要があります。

しかし単純に周辺化すると、\({\bf S}=\{s_1,s_2,...,s_N\}\)は観測データ数と同数あるため、計算量が爆発してしまいます。

この時、フォワード・バックワードアルゴリズムという手法が便利です。HMMのグラフィカルモデルの特徴の恩恵(\(x_t\perp s_{t+1}|s_t\), \(s_{t+1}\perp s_{t-1}|s_t\))を受け、計算量を削減することができます。ベイズ推論でお馴染みの「事後分布→結合分布→着目変数項以外を定数化」の流れで以下導出。\(f(s_n)\)と\(b(s_n)\)に分解できます。

$$q(s_n)=\sum_{{\bf S_{\backslash n}}}q({\bf S}) \propto  \sum_{{\bf S_{\backslash n}}}\tilde{p}({\bf X_{1:n}}, {\bf S_{1:n}})=f(s_n)b(s_n)$$$$f(s_n)=\tilde{p}(x_n|s_n)\sum_{{\bf S_{1:n-1}}}\tilde{p}(s_n|s_{n-1})\tilde{p}({\bf X_{1:n-1}},{\bf S_{1:n-1}})$$$$b(s_n)=\sum_{{\bf S_{n+1:N}}}\tilde{p}({\bf X_{n+1:N}},{\bf S_{n+1:N}}|s_n)$$

さらに詳しく分解すると、以下のように再帰式が得られます。計算量が\(O(N)\)に削減できています。 

  • フィードフォワード(前向き)$$f(s_n)=\tilde{p}(x_n|s_n)\sum_{{\bf S_{n-1}}}\tilde{p}(s_n|s_{n-1})f(s_{n-1})$$
  • バックワード(後向き)$$b(s_n)=\sum_{{\bf S_{n+1}}}\tilde{p}(x_{n+1}|s_{n+1})\tilde{p}(s_{n+1}|s_n)b(s_{n+1})$$ 

実装

Pythonにて実装しました。主にnumpyで書いています。データは自分で作成しました。

今回のコードをgithubに載せています。 githubはこちら 

 HMMの構造化変分推論のみ貼ります。

import numpy as np
import math

from scipy import special as sp
from scipy import stats as st


def HMM_paredict_s(x, K=3, iter_num=20, my_seed=0):
  #set seed
  np.random.seed(seed=my_seed)

  # sample num
  N=len(x)

  ### prior parameters
  # gamma distribution
  a=200
  b=5

  # Dirichlet distribution for \pi
  alpha = np.array([30,20,10])
  # Dirichlet distribution for \A
  beta = np.array([[20,10,10],[10,20,10],[10,10,20]])

### posterior parameters # gammma distribution a_update = np.array([200,200,200]) b_update = np.array([5,5,5]) # Dirichlet distribution for \pi alpha_update = np.array([30,20,10])
# Dirichlet distribution for \A beta_update = np.array([[20,10,10],[10,20,10],[10,10,20]]) # set state value first value s_mean=[] for n in range(N): s_mean.append([0.4,0.3,0.3]) s_mean = np.array(s_mean) for i in range(iter_num): ##################################################### # expectation of λ、π、A lam_mean=np.zeros(K) ln_lam_mean=np.zeros(K) ln_pi_mean=np.zeros(K) ln_A_mean=(np.zeros((K,K))).tolist() lam_mean = a_update/b_update ln_lam_mean=sp.digamma(a_update)-np.log(b_update) ln_pi_mean=sp.digamma(alpha_update)-sp.digamma(np.sum(alpha_update)) ln_A_mean = sp.digamma(beta_update) - sp.digamma(np.sum(beta_update, axis=0)) ##################################################### # q(Sn)、q(Sn, Sn-1) p_xn_sn_array =np.zeros((N, 3)) # forward algorithm f_pass = np.zeros((N, 3)) for n in range(N): ln_xn_sn = x[n]*ln_lam_mean-lam_mean p_xn_sn = np.exp(ln_xn_sn) p_xn_sn/=np.sum(p_xn_sn) p_xn_sn_array[n] = p_xn_sn p_s1 = np.exp(ln_pi_mean) p_s1/=np.sum(p_s1) p_sn_sn_1 = np.exp(ln_A_mean) p_sn_sn_1 /=np.sum(p_sn_sn_1, axis=0) if n==0: f_pass[n] = p_xn_sn * p_s1 else: f_pass[n] = p_xn_sn * np.sum(f_pass[n-1] * p_sn_sn_1, axis=1) f_pass[n] /=sum(f_pass[n]) # backward algorithm b_pass = np.zeros((N, 3)) for n in reversed(range(N)): if n==N-1: b_pass[n] = np.ones(3) else: ln_xn_sn = x[n+1]*ln_lam_mean-lam_mean p_xn_sn = np.exp(ln_xn_sn) p_xn_sn/=np.sum(p_xn_sn) p_sn_sn_1 = np.exp(ln_A_mean) p_sn_sn_1 /=np.sum(p_sn_sn_1, axis=0) b_pass[n] = np.sum(b_pass[n+1].reshape(3,1) * p_sn_sn_1 * p_xn_sn.reshape(3,1), axis=0) b_pass[n] /=sum(b_pass[n]) s_mean = f_pass * b_pass s_mean /=np.sum(s_mean, axis=1).reshape(N, 1) f_sn_1_b_sn = f_pass[1:].reshape(len(f_pass[1:]), 1, 3) * b_pass[:-1].reshape(len(b_pass[:-1]), 3, 1) s_s_1_mean = f_sn_1_b_sn * p_sn_sn_1 * p_xn_sn_array[1:].reshape(len(p_xn_sn_array[1:]), 3,1) s_s_1_mean = s_s_1_mean/np.sum(s_s_1_mean, axis=1).reshape(N-1,3,1) ########################################### # update a, b a_update = np.sum(x.reshape(len(x),1) * s_mean, axis=0) + a b_update = np.sum(s_mean, axis=0) + b # update α alpha_update = s_mean[0] + alpha # update β beta_update = np.sum(s_s_1_mean + beta, axis=0) ##################################################### # determine group number(just, symbol) by order of λ s_order = st.gamma(a=a_update, scale=1/b_update).mean().argsort() s_mean_ordered = s_mean[:,s_order] return s_mean_ordered

実験&結果

観測変数とモデルのみで潜在変数を予測しました。モデルで仮定した生成過程と実際のデータの生成過程が同じであるため当たり前ですが、潜在変数の正答率は94.5%でした。1列目が観測データ、2列目が正解の潜在変数、3列目が予測確率が最大になるカテゴリー、4列目が各カテゴリーの予測確率を表しています。

f:id:gashin_learning:20190828155130p:plain

パターン別のHMMとPMM(ポアソン混合モデル)比較実験

次に、潜在変数のマルコフ性を仮定しないポアソン混合モデル(PMM)との比較実験を紹介します。以下、パラメーターにも事前分布をもたしています。

f:id:gashin_learning:20190828164246p:plain

観測モデルである3つのポアソン分布の\(\lambda\)の値が近づいて差異を少なくしていった時のHMMとPMMの精度を比較しました。ポアソン混合モデルは完全分解変分ベイズで推論しています。

f:id:gashin_learning:20190828145109p:plain

下が結果です。

f:id:gashin_learning:20190828145132p:plain

図の通り、ポアソン分布の\(\lambda\)の値が近づいてくると、時系列情報のおかげでHMMの方が、精度が高く推論できています。
 

今回のコードをgithubに載せています。 githubはこちら

 

Twitterフォローよろしくお願いいたします! twitterはこちら