November 17, 2012

Pythonでパーティクルフィルタを実装してみる

パーティクルフィルタ(Particle filter)は,隠れマルコフモデルカルマンフィルタと同じように,システムの観測$Y$から状態$X$を推定する手法.どれもベイジアンフィルタに基づくもので,確率分布$p(X_t;Y_{0:t})$の表し方が異なる1のですが,パーティクルフィルタでは有限個のサンプリングによって確率分布を近似します.今回は重点サンプリング2を使ったパーティクルフィルタを実装してみます.ほかのフィルタに比べてループぐるぐる回すだけだからすごく簡単!


1. 隠れマルコフモデルはヒストグラム(離散),カルマンフィルタはガウシアン(パラメトリック),パーティクルフィルタはサンプリング(ノンパラメトリック)で表す
2. SciPyには有名ドコロの確率分布からサンプリングする関数が用意されている.任意の確率分布からサンプリングしたい場合には逆関数法,棄却サンプリング,重点サンプリングといった手法などを使う


パーティクルフィルタ


  • たくさんの粒子をばらまいておいて,それっぽく動かして,観測して,各々実際の観測とのズレを測って,正解っぽいっぽい粒子だけ残す,っていうのを繰り返す
  • 入力はN個のパーティクルの状態${\bf x}_t^{(i)}$,重み${\bf w}_t^{(i)}$ $(i=1,...,N)$と制御入力${\bf u}_t$と観測${\bf y}_t$
  • 出力は更新された状態${\bf x}_{t+1}^{(i)}$,重み${\bf w}_{t+1}^{(i)}$
  • 状態方程式$f$と観測方程式$g$が与えられている
    ${\bf x}_{t+1} = f({\bf x}_t, {\bf u}_t) + {\bf w} \leftrightarrow p({\bf x}_{t+1}|{\bf x}_t, {\bf u}_t)\\
    {\bf y}_t = g({\bf x}_t) + {\bf v} \leftrightarrow p({\bf y}_{t}|{\bf x}_t)$
  • 確率分布$p({\bf x}_t|{\bf y}_{0:t})$をN個の重み付きサンプル$\{w_t^{(i)}, {\bf x}_t^{(i)}\}$$(i=1,...,N)$で近似.$\delta$はデルタ関数.
    $p({\bf x}_{t}|{\bf y}_{0:t}) \approx \sum_{i=1}^{N} w_t^{(i)} \cdot \delta ({\bf x}_t - {\bf x}_t^{(i)})$
  • 新しい観測${\bf y}_{t+1}$があったら状態推定分布$p({\bf x}_{t+1}|{\bf y}_{0:t+1})$を3つのステップで更新

    1. 推定
    2. 更新
    3. リサンプリング


例題


2次元座標において、あるロボットが$t=0$に原点を出発して、速度$(4,4)$で動くとする。ロボットの進路は風などの影響を受け($\sigma_x=\sigma_y=2$),毎秒ごと4つの点$(0,0),(10,0),(0,10),(10,10)$からの距離を計測できて、計測には距離によらない誤差がある($\sigma_x=\sigma_y=4$)とする.このとき観測された軌跡から実際の軌跡を推定する.

Fig. 1 ピヨピヨ



推定


for i in range(N): ${\bf x}_{t+1}^{(i)} \sim p({\bf x}_{t+1}^{(i)}|{\bf x}_{t}^{(i)}, {\bf u}_{t})$
実際にはN個のパーティクルの位置を状態方程式に代入.
${\bf x}_{t+1}^{(i)} = f({\bf x}_{t}^{(i)}, {\bf u}_{t}) = {\bf A}{\bf x}_{t}^{(i)} + {\bf B}{\bf u}_{t} + {\bf w}$
ただし
${\bf A} = \left[
\begin{array}{cc}
1 & 0 \\
0 & 1 \\
\end{array}
\right],
{\bf B} = \left[
\begin{array}{cc}
1 & 0 \\
0 & 1 \\
\end{array}
\right],
{\bf w} \sim N(0, 2{\bf I})$


更新


for i in range(N): $w_{t+1}^{(i)} \leftarrow w_{t}^{(i)} \cdot p({\bf y}_{t+1}^{(i)}|{\bf x}_{t+1}^{(i)})$
尤度関数によって重みを更新.$\sum^i w_{t+1}^{(i)} = 1$で正規化.今回はモデルを正規分布にしたのでRBFカーネルで.尤度関数は推定値と観測値が似てれば似てるほど大きくなるように設定.物体追跡検知とかだと色の情報を使う.
$p({\bf y}_{t+1}^{(i)}|{\bf x}_{t+1}^{(i)}) \propto \exp(-\frac{(y-g(x))^2}{\sigma^2})$
ただし
$g({\bf x}) = \left[ ||{\bf x}-{\bf p}_1||, ||{\bf x}-{\bf p}_2||, ||{\bf x}-{\bf p}_3||, ||{\bf x}-{\bf p}_4|| \right]^{\mathsf{T}}\\
{\bf p}_1=\left[0, 0\right]^{\mathsf{T}}, {\bf p}_2=\left[10, 0\right]^{\mathsf{T}}, {\bf p}_3=\left[0, 10\right]^{\mathsf{T}}, {\bf p}_4=\left[10, 10\right]^{\mathsf{T}} \\
\sigma^2 = 4$


リサンプリング


$\{ {\bf x}_{t+1}^{(i)}, w_{t+1}^{(i)}=1 \} \leftarrow resampling(\{ {\bf x}_{t+1}^{(i)}, w_{t+1}^{(i)} \})$
重みに応じてリサンプリング.重みが偏らないように.毎回やる必要はない.色々手法があるらしいけど今回は単純に多項分布でサンプリング.


結果



Fig. 2 10ステップ分.が実際の軌跡.がパーティクル.がパーティクル平均.線形システムだから全然パーティクルフィルタ使う意味なかったけど…




実装


要NumPy

November 8, 2012

LDAの実装を試してみる

Latent Dirichlet allocationの実装を色々試してみた.自分でも実装したことある気がするけど.比較はまた後でやるとして使い方だけメモ.詳細は各リンク先で…
  1. Latent Dirichlet Allocation in C
  2. GibbsLDA++ A C C++ Implementation of Latent Dirichlet Allocation (LDA) using Gibbs Sampling for Parameter Estimation and Inference
  3. plda - A parallel C++ implementation of fast Gibbs sampling of Latent Dirichlet Allocation - Google Project Hosting

1. Latent Dirichlet Allocation in C


http://www.cs.princeton.edu/~blei/lda-c/

準備


lda-c-dist.tgz,ap.tgzをダウンロードしたら
tar xvfz lda-c-dist.tgz
tar xvfz ap.tgz
cd lda-c-dist
make

データ


各ドキュメントごとに行区切りで,語の種類数,語のインデックスと頻度を記述.インデックスはstringじゃないことに注意.
[M] [term_1]:[count] [term_2]:[count] ... [term_N]:[count]
たとえば
% head -n 3 ../ap/ap.dat
186 0:1 6144:1 3586:2 3:1 4:1 ...
174 68:1 512:1 514:2 3:1 4:1 ...
161 0:9 68:1 1538:1 3588:1 517:1 ...
1つ目のドキュメントは,語が186種類,語0が1回,語6144が1回…
結果を表示するためには番号と語を紐付けておくためのファイル(ap/vocab.txtみたいなやつ)も用意しておく.

実行と結果


- トピックの推定 estimation

以下を実行.結構時間かかるかも.
./lda est 1.0 50 settings.txt ../ap/ap.dat random test
testフォルダ以下に結果が出力されます.引数は,LDAのパイパーパラメータ$\alpha$,トピック数$K$,設定ファイル,データセット,初期状態,出力先.

$\alpha$に関しては $50/K$ にしておくといいらしい.$\beta$に関しては総語数に対する語彙数(異なり数)が多い場合は小さくするといいらしい.とどこかに書いてあった,気がする.ここでは$\beta$は0.1で固定っぽい.

- 他のデータセットの推定 inference

トピックの推定で使ったのと同じフォーマットのデータからディリクレパラメータと尤度を推定できる.意味ないけど上でトピックの推定をしたデータでディリクレパラメータの推定をする場合は以下.
./lda inf settings.txt test/final ../ap/ap.dat inference
出力はinference-gamma.dat,inference-lda-lhood.dat.

- 結果の表示

各トピックの上位10語を表示します.
python topics.py test/final.beta ../ap/vocab.txt 10



2. GibbsLDA++ A C C++ Implementation of Latent Dirichlet Allocation (LDA) using Gibbs Sampling for Parameter Estimation and Inference


http://gibbslda.sourceforge.net/

準備


GibbsLDA++: A C/C++ Gibbs Sampling LDA | Free Science & Engineering software downloads at SourceForge.netからGibbsLDA++-0.2.tar.gzをダウンロードして以下.
tar xvfz GibbsLDA++-0.2.tar.gz
cd GibbsLDA++-0.2
make clean
make all

データ


最初にドキュメント数,あとは行がドキュメントを表し,スペース区切りで語を羅列.
[M]
[document1]
[document2]
...
[documentM]
[documenti] = [wordi1] [wordi2] ... [wordiNi]
たとえば
% head -n3 trndocs.dat
1000
abil absenc acquisit acquisit agreem ...
activ ball ball band brief ...

実行と結果


Usageにある通り.

- パラメータ推定 estimation
src/lda -est -alpha 0.5 -beta 0.1 -ntopics 100 -niters 1000 -savestep 100 -twords 20 -dfile models/casestudy/trndocs.dat
estでLDAのパラメータを推定します.LDAのハイパーパラメータalpha,beta,トピック数ntopics,繰り返し回数niters,ステップsavestep,出力語数twords,データdfile.twordsを指定すると,各トピックの特徴語が出力されます.

- 途中のモデルから
src/lda -estc -dir models/casestudy/ -model model-01000 -niters 800 -savestep 100 -twords 30
estcで指定したモデルからパラメータを推定します.

- 他のデータの推定 inference
src/lda -inf -dir models/casestudy/ -model model-01800 -niters 30 -twords 20 -dfile newdocs.dat
infで作ったモデルから他のデータセットの推定をします.



3. plda - A parallel C++ implementation of fast Gibbs sampling of Latent Dirichlet Allocation - Google Project Hosting


http://code.google.com/p/plda/

準備

tar xvfz plda-3.0.tar.gz
cd plda
make lda infer

データ


行ごとにドキュメントを表し,語 頻度を繰り返す.
[word1] [word1_count] [word2] [word2_count] [word3] [word3_count] ...
たとえば
% head -n3 testdata/test_data.txt
concept 1 consider 1 global 1 entropy 1 go 1 ...
externally 1 global 1 dynamic 1 resistance 1 illustrated 1 ...
consider 1 chain 1 global 1 leads 1 go 1 ...

実行と結果


- 訓練

./lda --num_topics 2 --alpha 0.1 --beta 0.01 --training_data_file testdata/test_data.txt --model_file testdata/lda_model.txt --burn_in_iterations 100 --total_iterations 150
パラメータは上のものと似たようなもん.testdata/lda_model.txtに出力.出力のそれぞれの行は語のトピックの分布を表す.たとえば
% head -n3 testdata/lda_model.txt
concept 179.3 2.7
consider 921.98 0.02
global 296.3 180.7

- 表示

訓練で得られた結果を見やすく表示する.
python view_model.py testdata/lda_model.txt

- 他のデータセットの推定

./infer --alpha 0.1 --beta 0.01 --inference_data_file testdata/test_data.txt --inference_result_file testdata/inference_result.txt --model_file testdata/lda_model.txt --total_iterations 15 --burn_in_iterations 10
alphaは訓練の時と同じものを使いましょう.

- パラレル

http://code.google.com/p/plda/wiki/PLDAManual


...


- 論文と実装の紐付けがまとまってるところないかあなあ
- 論文書いている人はみんな実装も公開してほしいなああああああああ