February 17, 2012

Pythonで線形動的システムのパラメータ推定を実装してみるもなんかうまくいってない

線形動的システム(線形力学系,Linear dynamical systems; LDS)のEMアルゴリズムによるパラメータ推定を行いました.が,うまくいってない気がします.カルマンフィルタなどで扱う状態方程式(外部入力なし)と観測方程式において,状態と観測が与えられたときに,係数を求めるのが目的です.Ghahramaniらの論文1の手法を実装します.

状態方程式と観測方程式が
$\begin{align}
{\bf x}_{t+1}&=A{\bf x}_t + {\bf w}_t\\
{\bf y}_{t}&=C{\bf x}_t + {\bf v}_t\\
{\bf w}_t &\sim N(0,Q)\\
{\bf v}_t &\sim N(0,R)\\
\end{align}$
というように表され,$\{{\bf y}_1,...,{\bf y}_T\}$が与えられたときに,$A,C,Q,R$と初期状態の平均ベクトル$\pi_1$と分散共分散行列$V_1$を推定します.

論文がとても分かりやすいので,実装に必要なところだけ抜粋します.


Eステップ

以下のように定義します.
$\begin{align}
{\bf x}_t^{\tau}&=E({\bf x}_t|\{{\bf y}_1^{\tau}\})\\
V_t^{\tau}&=Var({\bf x}_t|\{{\bf y}_1^{\tau}\})\\
\end{align}$

Kalman filter
$\begin{align}
{\bf x}_1^0 &= {\bf \pi}_1\\
V_1^0 &= V_1\\
{\bf x}_t^{t-1} &= A {\bf x}_{t-1}^{t-1}\\
V_t^{t-1} &= A V_{t-1}^{t-1} A' + Q\\
K_t &= V_t^{t-1}C'(CV_t^{t-1}C'+R)^{-1}\\
{\bf x}_t^t &= {\bf x}_t^{t-1} + K_t({\bf y}_t-C{\bf x}_t^{t-1})\\
V_t^t &= V_t^{t-1} - K_t C V_t^{t-1}\\
\end{align}$

Kalman smoother
$\begin{align}
J_{t-1} &= V_{t-1}^{t-1}A'(V_t^{t-1})^{-1}\\
{\bf x}_{t-1}^T &= {\bf x}_{t-1}^{t-1} + J_{t-1}({\bf x}_t^T - A{\bf x}_{t-1}^{t-1})\\
V_{t-1}^T &= V_{t-1}^{t-1} + J_{t-1}(V_t^T - V_t^{t-1})J_{t-1}'\\
V_{T,T-1}^T &= (I - K_T C)AV_{T-1}^{T-1}\\
V_{t-1,t-2}^T &= V_{t-1}^{t-1}J_{t-2}' + J_t-1(V_{t,t-1}^T - AV_{t-1}^{t-1})J_{t-2}'\\
\end{align}$

これらを用いて,Mステップへの入力を求めます.
$\begin{align}
{\hat {\bf x}}_t & = {\bf x}_t^T \\
P_t & = V_t^T + {\bf x}_t^T {{\bf x}_t^T}' \\
P_{t,t-1} & = V_{t,t-1}^T + {\bf x}_t^T{{\bf x}_{t-1}^T}' \\
\end{align}$


Mステップ
$\begin{align}
C^{new} & = (\sum_{t=1}^T {\bf y}_t {\hat{\bf x}}_t')(\sum_{t=1}^T P_t)^{-1} \\
R^{new} & = \frac{1}{T}\sum_{t=1}^T({\bf y}_t{\bf y}_t' - C^{new}{\hat{\bf x}}_t{\bf y}_t') \\
A^{new} & = (\sum_{t=2}^T P_{t,t-1})(\sum_{t=2}^T P_{t-1})^{-1} \\
Q^{new} & = \frac{1}{T-1}(\sum_{t=2}^T P_{t} - A^{new}\sum_{t=2}^T P_{t-1,t}) \\
{\bf \pi}_1^{new} & = {\hat {\bf x}}_1 \\
V_1^{new} & = P_1 - {\hat {\bf x}}_1{\hat {\bf x}}_1'
\end{align}$
なおN個の観測列があるときは
$\begin{align}
{\bar {\hat {\bf x}}}_t & = \frac{1}{N}\sum_{i=1}^N {\hat {\bf x}}_t^{(i)} \\
V_1^{new} & = P_1 - {\bar {\hat {\bf x}}}_1{\bar {\hat {\bf x}}}_1' + \frac{1}{N}\sum_{i=1}^N [{\bar {\hat {\bf x}}}_1^{(i)}-{\bar {\hat {\bf x}}}_1][{\bar {\hat {\bf x}}}_1^{(i)}-{\bar {\hat {\bf x}}}_1]' \\
\end{align}$



実験

Wikipediaのカルマンフィルタの例にあるシステムで実験してみました.トロッコが摩擦なしの無限レール上でランダムな力によって動くシステムです.

データを生成するために用意したパラメータは
A
[[ 1.   0.1]
 [ 0.   1. ]]
C
[[ 1.  0.]]
Q
[[  2.50000000e-05   5.00000000e-04]
 [  5.00000000e-04   1.00000000e-02]]
R
[[ 1.]]
でした.

が,下のコードによって推定された結果は
A
[[ 0.62656567  0.64773696]
 [ 0.56148647  0.0486408 ]]
C
[[ 0.19605285  1.14560846]]
Q
[[ 0.15479875 -0.12660499]
 [-0.12665761  0.31297647]]
R
[[ 0.59216228]]
pi_1
[[-0.39133729]
 [-0.89786661]]
V_1
[[ 0.18260503 -0.07688508]
 [-0.07691364  0.20206351]]
でした.下のコードでは生成されるデータは毎回異なるので結果はその都度変わります.がどうもうまくいかない.


赤が真の位置,マゼンダが推定された位置.青が真の速度,シアンが推定された速度.


どこがおかしいのだろう(:D)| ̄|_ =3


1. Z. Ghahramani and G.E. Hinton. Parameter estimation for linear dynamical systems. University of Toronto technical report CRG-TR-96-2, 6, 1996.

February 1, 2012

Xcode 4でPythonの開発環境を整える


XcodeでPythonを開発することができるようです.はじめは少々面倒かもしれませんが.でもXcodeのPythonシンタックスハイライトはきれいな気もするし,Emacs風のキーバインドも使えるし,アリかも.

Python in Xcode 4 - Stack Overflow訳しただけ参考にしました.

  1. Xcode 4を開く.
  2. "File" → "New" → "New Project…"で新しいプロジェクトを作成.


  3. "Mac OS X"の"Other"を選択.
  4. "External Build System"を選択し,"Next"をクリック.


  5. プロダクト名と識別名(Company Identifier)を入力.
  6. "Build Tool"には,/usr/bin/pythonを入力して"Next"をクリック.ほかの場所にPythonをインストールしてたらそれを入力.
  7. どこに保存するか決めて"Create".


  8. メニューバーから"File" → "New" → "New File…"を選択.


  9. "Mac OS X"の"Other"を選択.
  10. "Empty"を選択し,"Next"をクリック.


  11. 拡張子".py"を含んだPythonファイルをプロジェクトのフォルダに保存し"Save"をクリック.


  12. メニューバーから"Product" → "Edit Scheme…"を選択.


  13. 左のコラムから"Run"を選択.
  14. "Info"タブから,"Executable"フィールドを選択し,"Other…"を選択.


  15. 手順6で選んだ場所に移動.⇧⌘G (shift + command + g)で移動する場所を指定できる.


  16. "Executable"から選択.
  17. "Debugger"フィールドは,"None"を指定.


  18. "Arguments"タブから"Base Expansions On"フィールドを選択,プロジェクトの名前を選択.
  19. "Arguments Passed On Launch"の"+"をクリック.
  20. $(SOURCE_ROOT)/実行したいファイル名 を記入.ファイルはプロジェクトフォルダの中にないと実行できない.もしほかの場所のものを実行したいときはフルパスを記入.スペースを含む場合はクォーテーションで囲む.


  21. "OK"をクリック.
  22. コーディングをはじめましょう!



  • ⌘Rで実行できます.
  • デフォルトでは出力は下のパネルに表示されます.
  • タブ幅,文字コードは右のパネルから変更可能です.
  • Pythonの自動インデントはできません.
  • 応用すれば他の言語の開発環境もつくれるかも.

Macbook Airの環境を整える覚書 その2

自分用.以前書いたMacBook Airの環境を整える覚書と合わせて.

Finder


サイドバーの文字サイズ変更


System Preferences > GeneralでSidebar icon sizeをSmallに.
ついでにRestore windows when quitting and re-opening appsはアンチェック.



Terminal


ログインシェルをzshに変更


System PreferencesからUsers & Groups
左下のロックを外して,Current Userを右クリック,Advanced Options...
Login shellを/bin/zshに.


行数の変更,終了後の動作の変更,metaキーの設定の変更


PreferencesからWindowでRowsを43に.



PreferencesからShellでWhen the shell exits:をClose if the shell exited cleanlyに.



PreferencesからKeyboardでUse option as meta keyにチェック.




キーボード


フルキーボードアクセスをすべてのコンポーネントで.ダイアログでタブとスペースを組み合わせて選択できるようになる.




IME


Input Sources


System Preferences > Language & Text > Input SourcesからDvorak - Qwerty ⌘,Kotoeri > HiraganaRomajiの3つを選択.ときどき気が向いたらDvorak←
ついでにKeyboard & Character Viewerにもチェック.
2つ以上のInput Sourceの切り替えは⌘長押しで.

Google 日本語入力


すばらしいんだけどDvorakとうまく共存できない…?



App Store


Xcode

開発環境.ダウンロード後,Install Xcodeを実行.

Stuffit Expander

解凍.

The Unarchiver

解凍.

Evernote

メモ.

Skitch

スクリーンキャプチャ.



パッケージ管理


Homebrew


MacPortsより楽な気がする.
/usr/bin/ruby -e "$(curl -fsSL https://raw.github.com/gist/323731)"
要Intel CPU,OS X 10.5 or higher,Xcode with X11,Java Developer Update.



フォント


Ricty


Incosolata + Migu 1M.Terminalに開発環境に.RictyをインストールするためにはFontForgeを入れる必要があります.
brew install fontforge
ln -s /usr/local/Cellar/fontforge/20110222/FontForge.app /Applications
2行目のシンボリックリンクの張り方はFontForgeインストール後に表示される.



テキストエディタ


CotEditor


Preferences > Syntax > Delay coloringにチェック



ウェブブラウザ


Chrome


Firefox

Opera

LaTeX

UpTeX.app

色々まとまって入ってる.

ghostscript

Homebrewインストール後
brew install ghostscript

TexShop

LaTeX書くときは便利.テキストエディタで書いちゃうことも多い.

ファイル共有

Dropbox

よほど大きくならない限り計算結果なんかもここに.持ってるすべての端末でzsh, vim, emacs, CotEditorの環境設定も共有.~/Dropbox/Settings以下に共有する環境設定を置いておく.下は環境共有をまとめて設定するシェルスクリプト.つったってシンボリックリンク貼ってるだけ.

Python

pip

Python用パッケージ管理.http://www.pip-installer.org/en/latest/index.html
easy_install pip
easy_installのかわり.アンインストールもできる.

ScipySuperpack


Lionの環境に合わせたscipyやら入れるために.