satomacoto: 2011

December 31, 2011

LaTeX形式で記述した数式の画像を埋め込むBloggerガジェット

Bloggerで数式を表示させるために以前つくったChrome Extensionを使っていたのですが，再利用性を考慮して，LaTeX形式で記述された数式をBloggerで表示するためのHTML/Javascriptガジェットを作ってみました．ガジェットの設定の仕方はこちらの記事の下の方に．jQueryが読み込まれている必要があります．

レンダリングには Infographics / Google Chart API の Mathematical Formulas を利用しています．

以下のコードを貼付けます．使い方はコード中のコメントの通りです．クラス名が他とバッティングするようならば，変更したほうが良いかもしれません．

ax^2 + bx + c = 0
x = \frac{-b \pm \sqrt{b^2 - 4ac}}{2a}

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]だけ持っている．小麦を育てるためには肥料F₁ [kg/km^2]と農薬P₁ [kg/km^2]が必要である．一方，大麦を育てるためには肥料F₂ [kg/km^2]と農薬P₂ [kg/km^2]が必要である．それぞれの売値が小麦S₁[yen/kg]，大麦S₂[yen/kg]だとすると，売上を最大化するためにはどのように農地を振り分ければいいでしょう？

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

Maximize:

S₁x₁ + S₂x₂

Subject to:

0 ≤ x₁ + x₂ ≤ L
0 ≤ F₁x₁ + F₂x₂ ≤ F
0 ≤ P₁x₁ + P₂x₂ ≤ P
x₁ ≥ 0, x₂ ≥ 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

実装

上の問題でS₁=3, S₂=2, L=5, F₁=1, F₂=3, F=10, P₁=2, P₂=1, P=9とした線形計画法を解いてみる．

Maximize:

3x₁ + 2x₂

Subject to:

0 ≤ x₁ + x₂ ≤ 5
0 ≤ x₁ + 3x₂ ≤ 10
0 ≤ 2x₁ + x₂ ≤ 9
x₁ ≥ 0, x₂ ≥ 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)

最適解はx₁=4.0 x₂=1.0のとき14であることがわかりました．

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

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

November 28, 2011

Gistで書いたコードを行番号つきでBloggerに埋め込む

Gistを使えば簡単なコード断片を管理して他の人と共有できるらしい（いまさら）．gitレポジトリなので色々gitっぽく使えるらしい．Gistで書いたコードをBloggerに埋め込んでみる．

Gistでコードを書いて保存して，

embedをクリックするとスクリプトのリンク表示されるのでコピーして，

HTMLモードで編集して貼り付ける．

すると以下のように表示されます．

デフォルトだと行番号は表示されないので．表示させるために https://gist.github.com/454771 を参考にBloggerのレイアウトに以下のHTML/Javascriptガジェットを追加しました．（※他の場所でjQueryが読み込まれている必要があります．~~なんか番号のふられ方微妙だから後で直すかな…~~直して行番号は選択されないようにしました）

追記

without jQuery

November 23, 2011

ヘッドホンのステレオミニプラグの交換

SENNHEISER HD228の右が聞こえなくなりました．そこで秋月電子でステレオ・ミニプラグを買ってきて修理しました．

やり方は「3.5ステレオミニプラグ製作例」とか「How to replace a headphone jack」とかで検索
ミニプラグは秋月で¥90
HD228はeBayで$36
リッツ線はライターで炙る
赤は右，青は左，金はGND

秋月電子 3.5mmΦステレオ・ミニプラグ（金属タイプ）
SENNHEISER 密閉型ヘッドフォン HD 228

before

after

November 14, 2011

無料でGrowlをソースからインストールする

Growl 1.3.1の通知がうるさいのでアップデートする．App Storeで$1.99で手に入れることができますが，ソースからコンパイルしてインストールしちゃう．1.3だけど．

Growl - Growl Mercurial Accessを参考に．必要な環境は以下．

OS X Lion - 10.7
Xcode
Mercurial

手順

Growl 1.2.2のアンインストール
アンインストーラは http://growl.info/downloads
ちなみに古いバージョンのインストーラもおなじところに．
Mercurial SCMのインストール http://mercurial.selenic.com/
Google Code repositoryからソースを取得 http://code.google.com/p/growl/source/checkout
ターミナルから

hg clone https://code.google.com/p/growl/
を実行すると，growlってディレクトリができてローカルコピーされます．止まってんじゃないかくらい時間かかった．
growlディレクトリからGrowl.xcodeprojを開く
ビルドターゲットをGrowl.appに
ビルドしようとするとエラーが起こるかも

hgRevision.hがないとかそれっぽいこととか言われても無視．

サインエラーに関してはGrowlのページには

A signing error. You will need to generate a self signed certificate, and can do so at Keychain Access.app -> application's name menu -> certificate assistant -> create a certificate. Be sure to use the same prefix as the certificate that you are receiving errors about.
ってあるけどめんどくさいので Code Signing Identity を Don't Code Sign に．2箇所ある．
ビルド（⌘B）
Growl.appを右クリックしてShow in FinderするとGrowl.appが見つかるのでアプリケーションフォルダに移動して実行．終わり．

November 10, 2011

Google Scholarの検索結果に+1ボタンを追加するChrome Extension

Google Scholarの検索結果にはGoogle +1ボタンがついていないのでボタンを追加するChrome拡張機能をつくった．はじめは自分用に文献をソーシャルブックマークするシステムをApp Engineでつくろうかと思ったけどめんどくさいし，gdd11jpで+1ボタンの簡単だよって聞いたので．

インストール
Chrome Web Store - +1 to Google Scholar

構成

manifest.json
content_script.js
jquery.js (http://jquery.com/からダウンロードしてリネーム）

- manifest.json

- content_script.js

参考にしたページ

Content Scripts - Google Chrome Extensions - Google Code
Content Scriptのつくりかた
plusone-button
+1ボタンのつけかた

備考

+1ボタンの読み込み遅いせいか知らんけど結果の表示が遅くてイライラする…たぶんこいつのせいでGoogle Readerももたついてんじゃないの…
だからデフォルトで検索結果に+1ボタンつけてほしいです
しかし拡張機能のアイコンとかタイルとかGoogleのもの勝手に使ってるの怒られるのかな

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*アルゴリズムの意味なし。）

August 16, 2011

SIMロック解除とケータイ補償お届けサービスあるいは外装交換

SIMロック解除をドコモに頼んでやってもらうと、ケータイ補償お届けサービス（1回目5,250円、2回目8,400円）あるいは外装交換（5,250円?）を利用して新しい携帯電話にしてもSIMロック解除をされた状態にして受け取ることができる、と言う事をDSの店員さんに確認してGalaxy S2のSIMロック解除をしました（3,150円）。設定画面等からSIMロック解除されているかどうか確認することはできない、らしい…

docomoでは「GALAXY S II SC-02C」、SAMSUNGでは「GALAXY S2(SC-02C)」という表記。どっちでもいいのか。

(追記) 無事海外（タイ、中国、マカオ）でも使えました。

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, u_t, y_t+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 ロボットの位置は(x₀,x₁)で、3点からの観測値(y₀,y₁,y₂)を得る。

状態方程式は線形、観測方程式は非線形となります。

$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

June 19, 2011

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

カルマンフィルタは、時間変化するシステムの、誤差のある離散的な観測から現在の状態を推定する手法。Wikipediaの記事（カルマンフィルター）がわかりやすい。

状態方程式と観測方程式が次のように与えられているとき

${\bf x}_{t+1}={\bf Ax}_{t}+{\bf Bu}_{t}+{\bf w}_{t+1}$ （状態方程式）
${\bf y}_{t}={\bf Cx}_{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})$ （フィルタ分布）

線形カルマンフィルタ(LKF; Linear Kalman Filter)は μ_t, Σ_t, u_t, y_t+1 を入力として、 μ_t+1, Σ_t+1を出力する。1ステップのプロセスは以下のとおり。

# prediction

${\bf \hat{\mu}}_{t+1} \leftarrow {\bf A\mu}_{t}+{\bf Bu}$ （現在の推定値）
${\bf \hat{\Sigma}}_{t+1} \leftarrow {\bf Q}+{\bf A\Sigma}_t {\bf A}^T$ （現在の誤差行列）

# update

${\bf \tilde{y}}_{t+1} \leftarrow {\bf y}_{t+1}-{\bf C\hat{\mu}}_{t+1}$ （観測残差）
${\bf S}_{t+1} \leftarrow {\bf C}{\bf \hat{\Sigma}}_{t+1}{\bf C}^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}{\bf \hat{\Sigma}}_{t+1}$ （更新された現在の誤差行列）

観測を得るごとにPredictionとUpdateを繰り返すことで、現在の状態を推定します。

導出は後述（予定）。

例題を。

2次元座標において、あるロボットがt=0に原点を出発して、速度(2,2)で動くとする。ロボットの進路は風などの影響を受け（σ_x=σ_y=1）、毎秒ごとに観測できるGPSによる位置座標には計測誤差（σ_x=σ_y=2）があるとする。このとき、観測された軌跡から実際の軌跡を推定する。

係数は

${\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 u}_{const}=\left[ \begin{array}{c} 2 \\ 2 \\ \end{array} \right]$ ${\bf Q}=\left[ \begin{array}{cc} 1 & 0 \\ 0 & 1 \\ \end{array} \right]$
${\bf C}=\left[ \begin{array}{cc} 1 & 0 \\ 0 & 1 \\ \end{array} \right]$ ${\bf R}=\left[ \begin{array}{cc} 2 & 0 \\ 0 & 2 \\ \end{array} \right]$