Showing posts with label Python. Show all posts
Showing posts with label Python. Show all posts

April 5, 2013

PythonでWebサイトのビデオを抽出してYouTube Data APIクライアントでプレイリストに登録する

YouTube系まとめサイトの動画を流しっぱなしにしておくために，以前herokuにRailsでpost-tube.satomacoto.comというのを自分用につくった．が，無料で使えるデータベースの制限を超えちゃってたので，WebサイトのYouTubeへのリンクを抽出してYouTubeのPlaylistに登録するPythonのスプリクトをやっつけで書いてみた．こんな風にプレイリストを作成する→http://www.youtube.com/user/stmct/videos．忘れたときのためにメモ．


Google Data API


Google Data APIを使ってYouTubeのプレイリスト作成や動画登録を行う．Pythonから操作するためにgdata-python-clientが用意されている．このサイトからダウンロードした gdata-...zip を解凍したら

sudo python setup.py install

でインストール．

使い方は

Developer's Guide: Python - YouTube — Google Developers


APIキーの取得


http://code.google.com/apis/youtube/dashboard/でAPIキーを取得する．


ClientLogin


ローカルで動かすのでClientLoginで認証する．認証はPlaylistの作成や動画の登録に必要．下記のスクリプトでは以下の部分を設定する．YouTubeのアカウント名は http://www.youtube.com/user/stmct だったらstmctにあたるところ．Googleを2段階認証にしている場合はアプリケーション固有のパスワードを生成する．

username = 'YouTubeのアカウント名'
email = 'APIキーを取得したGoogleのアカウント名'
password = 'Googleのパスワード'
developer_key = '取得したAPIキー'


コード


使い方

python youtube_gdata_create_playlist.py http://...

汚いコードだな…




...


  • NAVERみたいに次のページがある場合に対応するか
  • 人がまとめたものを使うってのはまずいのかな
  • post-tube.satomacoto.comってどうやってつくったんだっけかな
  • 連続再生だけが目的じゃなくてサイトごとにどんなビデオが貼り付けられているかとか自分の再生履歴やらお気に入りやらで何かしようとか思ってたんだっけ
  • プレイリストの履歴は取れないからな…
At April 05, 2013 Labels 0 Comments

March 31, 2013

Pythonの辞書をvalue値でソートするコードの実行時間の比較

get，lambda，itemgetter，zipを使ったvalue値でのソートを比較してみた．
以下コード．



timeitで計測．

>>> import timeit
>>> timeit.timeit(stmt='get.sort_test()', setup='import get', number=10000)
0.6642911434173584
>>> timeit.timeit(stmt='itemgetter.sort_test()', setup='import itemgetter', number=10000)
0.6961650848388672
>>> timeit.timeit(stmt='zip.sort_test()', setup='import zip', number=10000)
0.6995840072631836
>>> timeit.timeit(stmt='lamb.sort_test()', setup='import lamb', number=10000)
0.7502970695495605

ただし
  • zipの場合(value, key)のリストになる
  • getの場合keyのリストになる（value値は返さない）
に注意．

よく見かけるのはlambdaを使うものだが，同じ結果を得たかったらoperator.itemgetterを使うほうが少し早い．上位だけ取ってきたいときはgetが早いかも．

他にもやり方があるだろうか．
At March 31, 2013 Labels 0 Comments

December 15, 2012

政党と政策の距離の可視化

政党間の類似度の可視化で使ったデータを標準化して特異値分解して可視化．双対尺度法とか数量化理論III類とかコレスポンデンス分析とか言われている．たしか．

図 ◆政策，■政党．政策とそれに賛成する政党が近くにある…はず

以下コード

At December 15, 2012 Labels 0 Comments

政党間の類似度の可視化

日本政治.comの投票マッチングから各政党の政策に関する質問に対する姿勢から，政党間の距離を計算し可視化してみた．


ただし以下に注意

  • 各質問の重みは考慮していない
  • 距離の定義を変えればまったく異なる見え方に

さらに政党と政策との距離も可視化してみた→政党と政策の距離の可視化


以下コード


  • やっつけなのでてきとう
  • あとで双対尺度法試してみる
  • 一つ前に書いた要約プログラムで連立政権について考えてみる
At December 15, 2012 Labels 0 Comments

December 10, 2012

PythonでDocument Summarization based on Data Reconstruction (AAAI 2012)の実装

Zhanying He, Zhejiang University, et al., Document Summarization based on Data Reconstructionを実装してみる。AAAI-12 Outstanding Paper Awards。


概要


  • 従来の要約は冗長性を最小化するようなメイントピックを含む文を抽出することによって実現。
  • 本手法はオリジナルのドキュメント全体を再現できるような文集合を抽出して再構成。そのために抽出した文集合を評価するために再構成関数 reconstruction function の提案。
    • 線形的再構成 linear reconstruction。文の線形的な組み合わせによってドキュメントの近似。貪欲法 greedy strategy で最適化。
    • 非負線形的再構成 nonnegative linear construction 。文の線形的な組み合わせを足し算で再構成。乗算型重み更新 multiplicative updating で最適化。
  • 提案するフレームワークをDSDR (Document Summarization based on Data Reconstruction)と呼ぶ。
  • DUC 2006とDUC 2007で実験。ランダム、Lead、LSA、ClusterHITS、SNMFと比較。


DSDR


要約文がなるべくドキュメント全体の内容を含むようにする。再構成誤差（reconstruction error）を小さくするようにする。
  • ドキュメントの各文について重み付き語頻度ベクトル weighted term-frequency vector で表現。ステミングをしておいたり、ストップワードを取り除いたりしておく。候補文集合 the candidate set 。
  • 再構成関数により候補文集合から選ばれた文集合の評価。
  • 再構成誤差を最小にするような最適な組み合わせを探索。

コンセプト


まずオリジナルドキュメントと要約の再構成誤差を
    $$L({\mathbf v}_i - f(X; {\mathbf a}_i))$$
と考える。ただし、候補文集合$V$、要約文集合$X$、全単語数$d$として、$V=[{\mathbf v}_1,...,{\mathbf v}_n]^T$ where  ${\mathbf v}_i \in R^d$, $X=[{\mathbf x}_1,...,{\mathbf x}_m]^T$, $n > m$とする。${\mathbf v}_i$は語頻度ベクトル、また$f(\cdot)$を再構成関数とする。このとき再構成誤差を最小とするように$X, A$を定めれば目的関数は
    $$\min_{X,A} \sum_{i=1}^n ||{\mathbf v}_i - f(X; {\mathbf a}_i)||^2$$
となる。これに対して2つの手法で解く。

本論文では再構成関数は
$$f_i(X;{\mathbf a}_i)=\sum_{i=1}^m{\mathbf x}_j a_{ij}$$
と表し、候補文集合が
$${\mathbf v}_i\approx\sum_{i=1}^m{\mathbf x}_j a_{ij}$$
と選ばれた文集合の線形結合で近似されるとする。

詳細は論文参照。

線形的再構成 linear reconstruction

$$\begin{aligned}\min_{X, A} & \sum_{i=1}^n||{\mathbf v}_i-X^TA||^2+\lambda||{\mathbf a}_i||^2\\s.t. & X \subset V, |X|=m \\& A = [{\mathbf a}_1, \cdots, {\mathbf a}_n]^T \in R^{n\times m}\end{aligned}$$

非負線形的 nonnegative linear reconstruction

$$\begin{aligned}\min_{{\mathbf a}_i, \beta} & \sum_{i=1}^n\{ ||{\mathbf v}_i -V^T {\mathbf a}_i||^2+\sum_{j=1}^n\frac{a_{ij}^2}{\beta_j} \} + \gamma||\beta||_1 \\s.t. & \beta_j \ge0, a_{ij} \ge 0, {\mathbf a}_i \in R^n\end{aligned}$$



実装


論文に疑似コードが載っているのでPythonで実装。要NumPy。

V = [[1,0,0,0],
     [0,1,0,0],
     [1,1,0,0],
     [0,0,1,0],
     [0,0,1,1]]
というドキュメント（行が文、列が語）に対して、
  • 線形的再構成では3番目と5番目の[1,1,0,0], [0,0,1,1]という文が選ばれた。
  • 非負線形的再構成ではそれぞれに[ 0.49301097 0.49301097 0.6996585 0.49301097 0.70211699]という重みがつき、同様に3番目と5番目の文が選ばれやすいことを示している。
それぞれトレードオフパラメータは尖りやすさ。これが小さいと過学習しがち。





実験


あとで。


メモ


  • 報知的要約。（cf. 指示的要約）
  • LSAと似てる気がする。制約ついてるのがちょっと違うのかな。
  • おもしろい。でも実装してみただけ。あとでちゃんと読む。
  • テストの書き方がわからない。
  • 訳し方がわからない。

At December 10, 2012 Labels , 0 Comments

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

At November 17, 2012 Labels , , 0 Comments

July 4, 2012

Pythonでラベル伝搬法を試してみる

ネットワークの構造を予測解析のタスクにはノード判別とリンク予測があります．ノード判別問題は，幾つかのノードについてクラスラベルが与えられているとき，残りのクラスラベルが未知のノードに対してクラスラベルを予測する問題です．

ノード判別手法の最も簡単なもののひとつとしてラベル伝搬法という手法が知られています．ラベル伝搬法のアルゴリズムのひとつを文献1の801ページに基づいて実装してみました．なおラベル伝搬法については文献2の11章にまとまってました．

  1. 鹿島, グラフとネットワークの構造データマイニング, 電子情報通信学会誌 93(9), 797-802, 2010.
  2. Chapelle, O. et al., Semi-supervised learning, MIT Press, 2006.


ラベル伝搬法は「ネットワーク上で隣り合ったノードは同じクラスに属する」と仮定して未知のノードにラベルを振る半教師あり学習の手法．ここでは +1 と -1 の2種類のクラス判別問題．

  • ネットワークの構造は${\bf W}$で表す．${\bf W}$の$i,j$成分は$i$番目のノードと$j$番目のノードにリンクがある(1)か，ない(0)か．
  • クラスラベルはベクトル${\bf y}$で表す．ふられてないときは0．
  • 予測値はベクトル${\bf f}$で表す．それぞれ[-1,1]の連続値．

で隣り合ったノードの予測値が互いに近くなるように決定するための目的関数は以下．

$\begin{align}
J({\bf f})&=\sum_{i=1}^l(y^{(i)}-f^{(i)})^2 + \lambda \sum_{i<j}y^{(i,j)}(f^{(i)}-f^{(j)})^2 \\
&=||{\bf y}-{\bf f}||_2^2+\lambda {\bf f}^T{\bf L}{\bf f}
\end{align}$

ただし${\bf L}\equiv {\bf D}-{\bf W}$で${\bf D}$は${\bf W}$の各行の和を対角成分に持つ行列．λは1項目と2項目のバランスを取る定数．1項目は正解に近づけ，2項目は隣合うのノードの予測値を近づけます．

で，この最小化問題の解が

$({\bf I}+\lambda {\bf L}){\bf f}={\bf y}$

の解として求められます．


以下，てけとーにノードクラスとリンクを決めて試してみました．要scipy．




以下，結果をD3.jsで可視化しました．

予測前

ブルーとグリーンが予めクラスを与えているノード．オレンジのノードのクラスを予測します．

予測後

λ=1の結果です．node1とnode3はブルー，node8とnode9はグリーンに分けられました．ノードの色の濃さ（白か黒か）でどっちに近いか示しています．




At July 04, 2012 Labels , 0 Comments

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.

At February 17, 2012 Labels , 1 Comments

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の自動インデントはできません．
  • 応用すれば他の言語の開発環境もつくれるかも．
At February 01, 2012 Labels , 0 Comments

January 13, 2012

PythonでMapReduceっぽいものを実装してみる

satomacoto: 青空文庫のルビのマイニングをやり直した

MapReduceっぽいものを触ってみたかったのでPythonで実装．参考にしたのは

Writing An Hadoop MapReduce Program In Python @ Michael G. Noll
http://www.michael-noll.com/tutorials/writing-an-hadoop-mapreduce-program-in-python/

あんま意味ないけどdefaultdictを使ったものと比較．

コード

以下mapper.pyとreducer.pyのコード．ほぼコピペ．




以下defaultdictのコード．



実行結果

https://github.com/downloads/satomacoto/Playground/count.zip

% head count.txt
 obscene house ナンバー・ナイン 1
A l'odeur du soleil sur les lavandes douces. もうむらさきにうれているげな 1
Abaisse'〕 アベッセ 4
Autant de pluie autant de tristesse, Paris qui m'oppresse! くさくさするほどあめがふる 1
Aux figuiers qui 〔mu^riront〕, au vent qui passera, みなみのくにではいちじくが 1
Belle-vue de Tombeau ベル・ビュウ・ド・トンボウ 2
Bonjour Monsieur ボンジュール・ムッシュウ 1
But this fold flow'ret climbs the hill この花こそは山にも攀ぢよ 1
Cafe' カフエ 1
Cafe' カツフエ 1
% wc count.txt
  260652  834471 6590672 count.txt


速度比較

% time cat ruby_rev.txt | python mapper.py | sort -k1,1 | python reducer.py > count_mapreduce.txt
cat ruby_rev.txt  0.00s user 0.07s system 0% cpu 59.480 total
python mapper.py  7.74s user 0.06s system 13% cpu 59.485 total
sort -k1,1  79.92s user 0.43s system 86% cpu 1:33.18 total
python reducer.py > count_mapreduce.txt  10.48s user 0.12s system 11% cpu 1:33.17 total

% time cat ruby_rev.txt | python count_defaultdict.py | sort -k1,1 > count_defaultdict.txt 
cat ruby_rev.txt  0.00s user 0.07s system 1% cpu 6.228 total
python count_defaultdict.py  6.90s user 0.12s system 98% cpu 7.131 total
sort -k1,1 > count_defaultdict.txt  6.21s user 0.06s system 46% cpu 13.438 total

Hadoopではやってない．つーか入れてない．やってみるか…
At January 13, 2012 Labels 0 Comments

青空文庫のルビのマイニングをやり直した

################

「ヶケ」を漢字として扱っていなかったため，下のファイルでも取り切れていないようです．（bunさんにご指摘頂きました．ありがとうございます．）

################

ついカッとなって青空文庫からルビをマイニングしてみたものの

上のエントリでマイニングしたデータを少しちゃんと見たら，恥ずかしながらちゃんとルビが取れていないところがあったのでやり直しました．ごめんなさい．今度はあってるといいけど．ついでに"ファイル名\t文字\tルビ\n"だと後々使いにくいそうなので，
作品ID\t人物ID\t文字\tルビ\n
って形式にしました．

https://github.com/downloads/satomacoto/Playground/ruby_rev.zip

UTF-8．凡例含む出てきた順．被りあり．

################

「々仝〆〇ヶ」と「ケ」を漢字扱いにしてルビを抜き出しました．「ヶ原 はら」のような余計なものまで入っているので後処理が必要．

https://github.com/downloads/satomacoto/Playground/ruby_rev3.zip

################


% head ruby_rev.txt
046658 001257 倦怠 けんたい
046658 001257 玩味 がんみ
046658 001257 倦怠 けんたい
046658 001257 窪地 くぼち
046658 001257 鶉 うずら
046658 001257 啄木鳥 きつつき
046658 001257 叩 たた
046658 001257 漲 みな
046658 001257 栗鼠射 りすう
046658 001257 胡桃 くるみ
% wc ruby_rev.txt
 2212232 9099124 66071554 ruby_rev.txt

どうも
  • 青空文庫では漢字扱いの「々」が抜けてた
  • 注釈つきの文字（「※［＃「火＋稲のつくり」、第4水準2-79-87］みたいなの）」が取れてなかった
っぽいので，ルビを抜き出す正規表現を
((※［[^］]+?］|[一-龠々])+?|(?<=｜)([^｜]+?))《([^》]+?)》
としました．

どうなんだろうまだ微妙なのかな…正規表現ェ…
At January 13, 2012 Labels 1 Comments

January 7, 2012

ついカッとなって青空文庫からルビをマイニングしてみたものの

################

下のzipは上手いことルビが取れてないので以下を参照．
satomacoto: 青空文庫のルビのマイニングをやり直した

################

ついカッとなってWikipediaからカッコ表現をマイニングしてみた - nokunoの日記

を参考にさせていただき青空文庫のルビを抜き出してみました．なーんも使い道考えてないけど．MLやらNLPやらIMEやらの研究でもサーベイしてみて色々ゴニョゴニョしてみたい．まーでもとりあえずなんか適当に処理して辞書に登録してみて近代の文学者のごとく認《したた》めるかな．

https://github.com/downloads/satomacoto/Playground/ruby.zip

データは
ファイル名\t語\tルビ\n
って具合に出てきた順にタブ区切りで出力しました．語のカウントとかしてません．被りあります．凡例の分も出力しちゃってます．

% head ruby.txt
1000_ruby_2956.zip 小焦 こじ
1000_ruby_2956.zip 尾 ぴき
1000_ruby_2956.zip 小焦 こじ
1000_ruby_2956.zip 強請 ねだ
1000_ruby_2956.zip 葛飾 かつしか
1000_ruby_2956.zip 堤 どて
1000_ruby_2956.zip 雪洞 ぼんぼり
1000_ruby_2956.zip 堰 せ
1000_ruby_2956.zip 御留川 おとめがわ
1000_ruby_2956.zip 殺生 せっしょう
% wc ruby.txt
 1847315 5750935 65440323 ruby.txt

ルビ表記は 青空文庫 - Wikipedia にあるルールに従って抜き出しました．
ルビの表記は ｜と《》によって表現される。ルビを《》で囲んだり｜でルビのかかる文字列を特定するのは、視覚障碍者読書支援協会（BBA）の原文入力ルールに合わせたものである。

青空｜文庫《ぶんこ》
とあれば、「ぶんこ」というルビが「文庫」についていることを示す。

本日は晴天《せいてん》なり。
のように、仮名と漢字の間に｜が入る場合は｜を省略することも出来る。

｜ブルースカイ《青空》
のように、仮名にルビを強制的に振る時に使用することもある。

ちなみに「漢字《かんじ》」というパターン以外には
1029_ruby_20617.zip 接吻 ベエゼ
1040_ruby_20510.zip 然 きん／＼ぜん
1182_ruby_20549.zip 盡 こと／″＼
1067_ruby_1929.zip Esteros de Patino エステロス・デ・パチニョ
1317_ruby_22263.zip 三たび魔女の呪詛に萎れ毒気に染みぬる ウイズ・ヘキッツ・バン・スライス・プラステッド・スライス・インフェクテッド
1490_ruby.zip 願掛 がんがけ［＃底本では「け」が脱落］
15971_ruby_28461.zip 並木道 ブリ※［＃濁点付き片仮名ワ、1-7-82］ール
なんてパターンがあります．ほかにもパターンあると思う．漢字側に注があるのもあった．あとでなんとかするか．


以下抽出した手順のメモ．Lion．要wget（←brew install wget）．

手順メモ
  1. 青空文庫リストCSVのダウンロード

    青空文庫のこのページの"「公開中　作家別作品一覧拡充版：全て(CSV形式、zip圧縮）」をダウンロード"からリストをダウンロードして解凍するとCSVファイルが．Shift_JISなのでUTF-8に変換しておく（別に変換しなくても大丈夫だけど）．適当にディレクトリ掘ってから．
    wget http://www.aozora.gr.jp/index_pages/list_person_all_extended.zip
    unzip list_person_all_extended.zip
    iconv -f Shift_JIS -t UTF-8 < list_person_all_extended.csv > list_person_all_extended_utf-8.csv
  2. テキストファイルURLの抽出

    ダウンロードするためにテキストファイルのURLの列だけ抜き出して保存．
    python urls.py list_person_all_extended_utf-8.csv > urls.txt
  3. テキストファイルのダウンロード

    "zip"フォルダをつくって，そこにテキストのzipファイルをwgetでダウンロード．
    mkdir zip
    wget -i urls.txt -P zip
  4. ルビの抽出

    上のルールに従いルビを抽出．エラー起こっちゃったときは無視してるので，すべてのテキストに対してルビを取り出せているわけじゃない．
    python ruby.py 'zip/*.zip' > ruby.txt
  5. （おまけ）ルビがひらがな（と'・'と'ー'）だけのもの


わかんないこと

  • そもそもどうやって辞書に登録するんだろう？まとめてできんのかな？品詞どーすればいい？つーか既に辞書にある語もあるような．全部登録していいのかな？記号関連どうすんだろ？
  • 文語体⇔口語体，何とかできんのかな？音便やらの仮名遣パターンあるのかな？
  • 普通は「漢字《かんじ》」だけどそれ以外のも結構ある．おもしれーからむしろこっちをゴニョゴニョしてみるかな？
  • …


IME関連の研究おもしろそうだよなー
At January 07, 2012 Labels 0 Comments

January 6, 2012

Pythonでカーネル主成分分析を実装してみる

カーネル主成分分析(Kernel principal component analysis; kernel PCA)はその名のとおりカーネルを使って主成分分析(PCA)を行う手法．カーネルっていうのは距離みたいなもん（→Kernel trick - Wikipedia）．線形主成分分析はデータが元の変数のベクトルの線形結合で表されてるとしているけど，そうでもないときもあるよねってことで，非線形写像して主成分分析する．分散共分散行列の代わりに，各データ間のカーネルを要素とするカーネル行列の固有値固有ベクトルを求めるってこと．

Kernel pricipal component analysis - Wikipediaに載ってる画像っぽいのをつくってやってみた．つくったデータ→
kpca.data
"x1 x2 class\n"．3つのクラス．各クラス100個ずつ．

元データの描画．


カーネルは2通り試してみた．固有値大きい方から2つ取り出して固有ベクトルの要素の値をプロット．

その1 多項式カーネル
$$k({\mathbf x},{\mathbf y}) = ({\mathbf x}'{\mathbf y} + 1)^2$$



その2 ガウシアンカーネル
$$k({\bf x},{\bf y}) = \exp \left( \frac{- ||{\bf x} - {\bf y}|| ^2}{2\sigma^2}\right)$$



以下実装．要numpy，matplotlib．



距離（非類似度）行列の固有値固有ベクトル求めて描画するのはMDSっぽい
At January 06, 2012 Labels , , 0 Comments

January 5, 2012

Pythonで多次元尺度構成法を実装してみる

多次元尺度構成法(Multidimensional scaling; MDS)は多変量解析の手法です．よくデータ間の（非）類似度の情報可視化に用いられます．MDSは基本的には似ているものは近くに，似ていないものは遠くに配置するような座標を求めます．ここでは古典的（計量的; metric）MDSをPythonで実装してみます．

古典的MDSは以下の手順で可視化．実装に必要なところのみ．

  • 要素の値が距離の2乗と見なせる非類似度行列Sを用意する
  • Sにヤング・ハウスホルダー変換を施してPとする
  • Pをスペクトル分解する（固有値・固有ベクトルを求める）
  • 固有値の大きい方から2~3個選び，対応する固有ベクトルを取り出す
  • 各固有ベクトルの要素値をプロットする（2個の固有ベクトルの時は2次元，3個の固有ベクトルの時は3次元）

一応，下部に参考リンクを挙げましたが，詳しい説明は検索すれば山ほど…．ヤング・ハウスホルダー変換は距離の2乗の行列に両側から中心化行列をかける演算のことらしいです．中心化行列は距離の基準の原点を重心に移動するための行列らしいです．

ヤング・ハウスホルダー変換:
SがN×N行列のとき
$$P = -\frac{1}{2}({\bf I}-\frac{1}{n}{\bf 11}')S({\bf I}-\frac{1}{n}{\bf 11}')$$
と計算します．1は要素がすべて1のN次ベクトル[1,...,1]'です．なので11'は要素がすべて1のN×N行列．

アメリカの空港間の距離が与えられているときMDSで可視化してみます．
atl   chi   den   hou    la    mi    ny    sf   sea    dc
atl     0
chi   587     0
den  1212   920     0
hou   701   940   879     0
 la  1936  1745   831  1374     0
 mi   604  1188  1726   968  2339     0
 ny   748   713  1631  1420  2451  1092     0
 sf  2139  1858   949  1645   347  2594  2571     0
sea  2182  1737  1021  1891   959  2734  2408   678     0
 dc   543   597  1494  1220  2300   923   205  2442  2329     0

以下コード．要numpy, matplotlib．



結果はこんなかんじ．ラベルは手作業で付けました．


ちなみに実際の位置はこんなかんじ．

View US 10 Airport in a larger map

結果を180度回転すればだいたいなんとなくあってるかな？

上の例では非類似度行列を実際の距離の行列としましたが，MDSはデータ間の非類似度を測ることができたらそれっぽく2次元(3次元)にマッピングできるよ，ってことなので色々使えるかも知れません．

ちなみにRだったらcmdscaleでMDSが実装されてて楽

At January 05, 2012 Labels , , 0 Comments

December 3, 2011

Pythonで線形計画法を解いてみる


線形計画法（Linear Programming; LP or linear optimization）は，目的関数も制約関数も線形であるときの最適化問題を解く手法です．SciPy + GLPK + CVXOPT + OpenOpt + FuncDesignerで実装されているのを使ってみました．説明はWikipediaから抜粋．

線形計画法は以下のようにある制約関数のもとである目的関数を最大化あるいは最小化するもの．

maximize or minimize （目的関数）
{\b c}^T {\b x}
subject to （制約関数）
{\mb A}{\b x} \le {\b b}
{\b x} \ge {\b 0}
xは変数，A,b,cは係数．制約関数の形が凸多面体であるので目的関数は最適化可能．

たとえば
ある農夫がL [km^2]の広さの農地を持っているとき小麦と大麦を育てて売上を最大化させたい．農夫は肥料をF [kg]と農薬をP [kg]だけ持っている．小麦を育てるためには肥料F1 [kg/km^2]と農薬P1 [kg/km^2]が必要である．一方，大麦を育てるためには肥料F2 [kg/km^2]と農薬P2 [kg/km^2]が必要である．それぞれの売値が小麦S1[yen/kg]，大麦S2[yen/kg]だとすると，売上を最大化するためにはどのように農地を振り分ければいいでしょう？

という問題を線形計画法で表すと

Maximize:
S1x1 + S2x2
Subject to:
0 ≤ x1 + x2 ≤ L
0 ≤ F1x1 + F2x2 ≤ F
0 ≤ P1x1 + P2x2 ≤ P
x1 ≥ 0, x2 ≥ 0

これを行列形式で書くと

Maximize
\left[\begin{array}{cc}S_1 & S_2 \\\end{array}\right]\left[\begin{array}{c}x_1 \\x_2\end{array}\right]
Subject to

\lef[\beg{a}{c}x_1\\x_2\end{a}\right]\ge\lef[\beg{a}{c}0\\0\end{a}\right]

となる．



準備


環境はMac OS X 10.7 (Lion)．要XCode 4.2，Git，Homebrew，setuptools．


SciPy http://fonnesbeck.github.com/ScipySuperpack/

以下を実行．
git clone git://github.com/fonnesbeck/ScipySuperpack
cd ScipySuperpack
sh install_superpack.sh
Are you installing from a repository from this machine?にはyと答える．


GLPK (GNU Linear Programming Kit) http://www.gnu.org/software/glpk/

GLPK(Gnu Linear Programming Kit)は，線形計画問題(LP)や混合整数計画問題(MIP)を解くためのソルバー．GNU GPL．Homebrewを使ってインストール．
sudo brew install glpk


CVXOPT http://abel.ee.ucla.edu/cvxopt/download/index.html

CVXOPTはPythonの凸最適化(convex optimization)のためのパッケージ．sageでも利用可．GNU GPL．左のメニューのDownloadを選び，Installation from sourceのPython 2.5+をダウンロードして./srcのsetup.pyのBUILD_GLPKを0から1に書き換え，
# Set to 1 if you are installing the glpk module.
BUILD_GLPK = 1
以下を実行．
sudo setup.py install


OpenOpt + FuncDesigner http://www.openopt.org/Install

OpenOptは数値最適化フレームワーク．自前のソルバーに加え，複数のソルバーが使える．BSD．以下を実行．
sudo easy_install openopt
sudo easy_install FuncDesigner



実装


上の問題でS1=3, S2=2, L=5, F1=1, F2=3, F=10, P1=2, P2=1, P=9とした線形計画法を解いてみる．

Maximize:
3x1 + 2x2
Subject to:
0 ≤ x1 + x2 ≤ 5
0 ≤ x1 + 3x2 ≤ 10
0 ≤ 2x1 + x2 ≤ 9
x1 ≥ 0, x2 ≥ 0

コード

結果
------------------------- OpenOpt 0.36 -------------------------
solver: glpk problem: unnamed type: LP goal: max
iter objFunVal log10(maxResidual)
0 -0.000e+00 -100.00
GLPK Simplex Optimizer, v4.47
3 rows, 2 columns, 6 non-zeros
Preprocessing...
3 rows, 2 columns, 6 non-zeros
Scaling...
A: min|aij| = 1.000e+00 max|aij| = 3.000e+00 ratio = 3.000e+00
Problem data seem to be well scaled
Constructing initial basis...
Size of triangular part = 3
* 0: obj = 0.000000000e+00 infeas = 0.000e+00 (0)
* 2: obj = -1.400000000e+01 infeas = 0.000e+00 (0)
OPTIMAL SOLUTION FOUND
1 1.400e+01 -100.00
istop: 1000 (optimal)
Solver: Time Elapsed = 0.01 CPU Time Elapsed = 0.006718
objFunValue: 14 (feasible, MaxResidual = 0)
(4.0, 1.0)
最適解はx1=4.0 x2=1.0のとき14であることがわかりました．

確認のためにGoogleを使ってグラフで描画してみます．制約条件が青，赤，黄．目的関数が緑．最大化なら目的関数はなるべく右上の方に持っていきたい．ちょうど青と黄の交点が最適解(4,1)になっていることがわかります．


今回はだいぶ簡単な例題を扱いましたが，OpenOptではもっと大規模なもの扱えるようです．またOpenOptでは線形計画法のみならず，二次計画法(QP)や非線形計画法(NLP)など様々な最適化問題が解けるようです．CookbookにOpenOptやFuncDesingerで扱えるサンプルがまとまっているようです．あ，っていうかSciKitのほうを使ってみればよかったか…？
At December 03, 2011 Labels , 4 Comments

August 23, 2011

PythonでA*アルゴリズムを実装してみる

PythonでNexworkXを使ってA*アルゴリズムを実装してみます。まあNetworkXではA*アルゴリズムは実装されているのですが。

NetworkXは複雑ネットワークでいろいろやるためのPythonのパッケージ。easy_installでインストールするのが楽。

A*（A-star, エースター） アルゴリズムは、経路探索アルゴリズムの一つ。Wikipediaの解説がわかりやすい。(→Wikipedia ja, Wikipedia en)

ここではWikipediaにある疑似コードをほぼそのままの形で実装してみます。ランダムな重みつきエッジで構成される3x3のグリッドにおいて(0,0)から(2,2)への最短ルートを探索します。

以下ではヒューリスティック関数h*(n)として直線距離を使います。
が、上の設定では直線距離は
\forall n,0\leq h^*(n) \leq h(n)
を満たさないので、最短経路であることは保証されません…。

なので、一応、ヒューリスティック関数h*(n)による推定値を常に0とすることで、ダイクストラ法による最短距離を求めてみます。（A*アルゴリズムの意味なし。）

At August 23, 2011 Labels 0 Comments

June 22, 2011

Pythonで拡張カルマンフィルタを実装してみる

拡張カルマンフィルタ(EKF; Extended Kalman Filter)は非線形カルマンフィルタのひとつ。線形カルマンフィルタは線形システムを対象としていましたが、拡張カルマンフィルタは非線形システムを対象とします。ナビゲーションやらGPSやらに利用されている、らしい。

システムが非線形のとき、すなわち
{\bf x}_{t+1}=f({\bf x}_{t},{\bf u}_{t})+{\bf w}_{t+1} （状態方程式）
{\bf y}_{t}=h({\bf x}_{t})+{\bf v}_{t} （観測方程式）
{\bf w}_{t}\sim N(0,{\bf Q}){\bf v}_{t}\sim N(0,{\bf R}) （ノイズは正規分布）
p({\bf x}_{t}|{\bf y}_{0:t})=N({\bf x}_{t};{\bf \mu}_{t},{\bf \Sigma}_{t}) （状態は正規分布）
とするとき、関数 f は前の状態から推定値を与え、関数 h は観測値を与えますが、どちらの関数も直接共分散を求めることはできません。が、拡張カルマンフィルタでは状態方程式も観測方程式も微分可能であれば線形である必要はありません。

拡張カルマンフィルタでは状態方程式と観測方程式の線形化をするために、線形カルマンフィルタにおける時間遷移モデルと観測モデルに各関数の偏微分行列（ヤコビアン）を用います。
{\bf A}_{t} = \left\ \frac{\partial f({\bf x}_{t},{\bf u}_{t})}{\partial {\bf x}_{t}} \right|_{{\bf \mu}_t, {\bf u}_t}
{\bf C}_{t} = \left\ \frac{\partial h({\bf x}_{t})}{\partial {\bf x}_{t}} \right|_{{\bf \mu}_t}
あとは、線形カルマンフィルタと同じ。 μt, Σt, ut, yt+1 を入力として、 μt+1, Σt+1を出力します。1ステップのプロセスは以下のとおり。

# prediction
{\bf \hat{\mu}}_{t+1} \leftarrow f({\bf \mu}_{t},{\bf u}) （現在の推定値）
{\bf A}_{t} \leftarrow \left\ \frac{\partial f({\bf x}_{t},{\bf u}_{t})}{\partial {\bf x}_{t}} \right|_{{\bf \mu}_t, {\bf u}_t} （状態方程式の現在のヤコビアン）
{\bf \hat{\Sigma}}_{t+1} \leftarrow {\bf Q}+{\bf A}_t{\Sigma}_t {\bf A}_t^T （現在の誤差行列）
# update
{\bf C}_{t+1} \leftarrow \left\ \frac{\partial h({\bf x}_{t})}{\partial {\bf x}_{t}} \right|_{{\bf \hat{\mu}}_{t+1}} （観測方程式のヤコビアン）
{\bf \tilde{y}}_{t+1} \leftarrow {\bf y}_{t+1}-{\bf C}_{t+1}{\hat{\mu}}_{t+1} （観測残差）
{\bf S}_{t+1} \leftarrow {\bf C}_{t+1}{\bf \hat{\Sigma}}_{t+1}{\bf C}_{t+1}^T+{\bf R} （観測残差の共分散）
{\bf K}_{t+1} \leftarrow {\bf \hat{\Sigma}}_{t+1}{\bf C}^T{\bf S}^{-1} （最適カルマンゲイン）
{\bf \mu}_{t+1} \leftarrow {\bf \hat{\mu}}_{t+1}+{\bf K}_{t+1}{\bf \tilde{y}}_{t+1} （更新された現在の推定値）
{\bf \Sigma}_{t+1} \leftarrow {\bf \hat{\Sigma}}_{t+1}-{\bf K}_{t+1}{\bf C}_{t+1}{\bf \hat{\Sigma}}_{t+1} （更新された現在の誤差行列）


例題。
2次元座標において、あるロボットがt=0に原点を出発して、速度(2,2)で動くとする。ロボットの進路は風などの影響を受け（σx=σy=1）、毎秒ごと3つの点(0,0),(10,0),(0,10)からの距離を計測できて、計測には距離によらない誤差がある（σx=σy=2）とする。このとき、観測された軌跡から実際の軌跡を推定する。
Fig.1 ロボットの位置は(x0,x1)で、3点からの観測値(y0,y1,y2)を得る。

状態方程式は線形、観測方程式は非線形となります。
f({\bf x},{\bf u})=\left[\begin{a}{cc}1&0\\0&1\end{}\right]{\bf x}+\left[\begin{a}{cc}1&0\\0&1\end{}\right]{\bf u},\ \ \ {\bf u}=\left[\begin{a}{c}2\\2\end{}\right]
h({\bf x})=\left[\begin{array}{c} \sqrt{x_0^2+x_1^2}\\ \sqrt{(x_0-10)^2+x_1^2}\\ \sqrt{x_0^2+(x_1-10)^2} \end{array}\right]

観測方程式のヤコビアンは
J=\left[\begin{array}{cc} \fr{\pa{h}_0}{\pa{x_0}}&\fr{\pa{h}_0}{\pa{x_1}} \\ \fr{\pa{h}_1}{\pa{x_0}}&\fr{\pa{h}_1}{\pa{x_1}} \\ \fr{\pa{h}_2}{\pa{x_0}}&\fr{\pa{h}_2}{\pa{x_1}} \\ \end{array}\right]
として求められます。実装では上式を解析的に解く場合、数値的に解く場合を試してみます。なおSciPyには数値的にヤコビアンを求める関数が用意されているので、ここではそれを使ってみます。（微小区間ズラして傾き求めている…）

以下、実装。61行目と62行目で解析的か数値的か切り替えます。



割とうまくいっている実行結果。

Fig. 2 赤が実際の軌跡。青が推定値。
左が解析的なヤコビアン、右が数値的なヤコビアンを使った結果。

なんかうまくいってないような(?)結果も出てしまいます。

Fig. 3 なんか微妙…

UKFも実装してみようかな。

http://docs.scipy.org/doc/scipy/reference/optimize.html
http://old.nabble.com/calculating-numerical-jacobian-td20506078.html
http://projects.scipy.org/scipy/browser/trunk/scipy/optimize/slsqp.py
http://en.wikipedia.org/wiki/Broyden's_method
At June 22, 2011 Labels , 2 Comments
Subscribe to: Comments (Atom)