行列指数関数の公式

行列指数関数の公式をいつも忘れるので,ここに書き留めておく.

ブロック行列
{ \displaystyle
A = \begin{bmatrix}
A_{11} & A_{12}\\0&A_{22}
\end{bmatrix}
}
に対して,
{ \displaystyle
e^{At} = \begin{bmatrix}
e^{A_{11}t} & \int_0^t  e^{A_{11}(t-\tau)}A_{12}e^{A_{22}\tau}\mathrm{d}\tau\\0&e^{A_{22}t}
\end{bmatrix}
}
が成り立つ.

この公式の重要なことは,右辺右上に行列指数関数の積分がでてくることである.
このかたちの積分はいろいろなところで現れる.たとえば,連続時間線形システム
{\displaystyle
\dot x = Ax+Bu
}
を入力をゼロ次ホールドして得られる離散時間システム
{\displaystyle
x_{k+1} = A_dx_k + B_du_k
}
について,
{\displaystyle
A_d = e^{AT}, \quad B_d = \int_0^T e^{A(T-\tau)}B\mathrm{d}\tau
}
とかける.ただし, Tはサンプリング周期である.これらの行列は,ブロック行列
{ \displaystyle
M = \begin{bmatrix}
A & B\\0&0
\end{bmatrix}
}
に対する行列指数関数 e^{MT}の部分行列(左上と右上)である.

Tex (TikZ) でブロック線図を書く

制御系の研究をしていると,ブロック線図を書く場面が多くある.
今回は,ブロック線図をTexのTikZを使って書いてみた.
単純なブロック線図であれば,PowerPointとかで書くよりも簡単な気がする.
線が若干斜めで美しくない,とかいうことが生じないのも気持ちがいい.
# TikZ初心者なので,もっと効率よく書けるとかがあれば教えていただければ幸いです.

本記事では,

\usepackage{tikz}
\usetikzlibrary{calc}
\usetikzlibrary{positioning}
\usetikzlibrary{arrows.meta}

を前提とする.

まず,ブロックや矢印などの設定.

[
% 線と文字の間の間隔を設定
every node/.style={outer sep=0.12cm, inner sep=0},
% 矢印の設定
arrow/.style={-{Stealth[length=0.25cm]}, thick},
% ブロック
block/.style={rectangle, draw, minimum height = 1cm, 
minimum width=1.2cm, thick, outer sep = 0},
% 加え合わせ点
sum/.style={thick, circle, draw, inner sep=0,
minimum size=6pt, outer sep=0},
% 引き出し点
point/.style={radius=2pt}
]

このうえで,

\node [block] (K){$K$};
\node [block, right=1 of K] (G){$G$};
\node [sum, left=1of K] (sum){};

とやれば,ブロック間を1cmあけてブロックが並ぶ.
f:id:kawa0616:20180326062943p:plain:w300


これらを矢印でつなぐには,

\draw[arrow] (sum) -- (K);

とすれば良い.信号の名前を書き込みたければ,

\draw[arrow] (K) -- (G) node [above, pos=0.5] {$u$};

という感じで書けばよい.aboveが信号線の上,pos=0.5は真ん中を表している.
ちなみに,pos=0.5はmidway に置き換えることもできる.

さらに,始点,終点がブロックでない矢印も,

\draw[arrow] (G.east) -- +(2, 0) node[right]{$y$};
\draw[arrow] (sum.west)+(-1, 0) node[left]{$r$} -- (sum.west)
node[below, xshift=-10pt]{$+$};

のような感じで書くことができる.
xshiftは横方向の調整で,これがないと矢印の終点の真下になってしまう.

f:id:kawa0616:20180326063600p:plain

最後に,フィードバックループを作る.
まずは,引き出し点.

\fill [point] (G.east)+(1, 0) circle coordinate (y);

coordinateで引き出し点の座標に名前を付けている.あとは,

\draw [arrow] (y) -- +(0, -1) -| (sum) node[left, yshift=-15pt] {$-$};

とすればフィードバックループを書くことができる.
"-|" で折れ線が書けるのが便利.
f:id:kawa0616:20180326064227p:plain

すべてをまとめると,次のようになる.

\begin{tikzpicture}
[
% 線と文字の間の間隔を設定
every node/.style={outer sep=0.12cm, inner sep=0},
% 矢印の設定
arrow/.style={-{Stealth[length=0.25cm]}, thick},
% ブロック
block/.style={rectangle, draw, minimum height = 1cm,
minimum width=1.2cm, thick, outer sep = 0},
% 加え合わせ点
sum/.style={thick, circle, draw, inner sep=0,
minimum size=6pt, outer sep=0},
% 引き出し点
point/.style={radius=2pt}
]
\node [block] (K){$K$};
\node [block, right=1 of K] (G){$G$};
\node [sum, left=1of K] (sum){};
\draw[arrow] (sum) -- (K);
\draw[arrow] (K) -- (G) node [above, pos=0.5] {$u$};
\draw[arrow] (G.east) -- +(2, 0) node[right]{$y$};
\draw[arrow] (sum.west)+(-1, 0) node[left]{$r$} -- (sum.west)
node[below, xshift=-10pt]{$+$};
\fill [point] (G.east)+(1, 0) circle coordinate (y);
\draw [arrow] (y) -- +(0, -1) -| (sum) node[left, yshift=-15pt] {$-$};
\end{tikzpicture}


頑張れば,もっと複雑なブロック線図もかける.
f:id:kawa0616:20180326064830p:plain
この例では,多入力のブロックがあったり,ちょっと複雑な配置も含まれている.
このブロック線図を出力するtexのソースは次の通り.

Bash on Windowsでシェルスクリプトを動かすときにハマった話

Windowsで作成したシェルスクリプトBash on Windowsで実行しようとして,ハマりまくった.

mkdirしようとして, 「Directory nonexistent」とか言われて意味不明だったが,
なんということはない,WindowsBash on Windowsで改行コードが異なることが原因だった.

ちなみに,改行コードはnkfコマンドでなおすことができる.
Windows => Bashならば,次の通り.

nkf -Lu dos.txt > dos2unix.txt

今回はBash on Windows上でこれを実行して,解決.

藤井四段のレーティングをベイズ推定してみた

藤井聡太四段が強すぎて,実際のレーティングがどれくらいなのか気になったのでベイズ的に計算してみた.

ベイズの定理から,

p(A|B) = \frac{p(A)p(B|A)}{\int p(A)p(B|A)\mathrm{d}A}

が成り立つ.
Aを藤井四段のレーティング,Bを2017年8月16日までの勝敗データとして,レーティングの分布p(A|B)を調べてみる.
これは,これまでの勝敗データから推定したときに,レーティングがどのあたりにありそうか,というのを確率密度で表現したものである.

ベイズの定理を用いるためには,事前分布 p(A)が必要である.
今回は簡単に0から5000の間の一様分布
 p(A) = \begin{cases}1/5000 & 0 \leq A \leq 5000\\ 0 & \mathrm{otherwise}\end{cases}
とする.すなわち,勝敗のデータがないときは,0から5000までのどこにあるかは全く予想できない,ということを仮定する.

また,尤度関数 p(B|A)については, i番目の対戦( i=1, \dots, N)に関する尤度関数
 p(B_i|A) = \begin{cases}\frac{1}{1+10^{(R_i-A)/400}}&(B_i=\mathrm{Win})\\\frac{1}{1+10^{(A-R_i)/400}}&(B_i=\mathrm{Lose})\end{cases}
を用いて,
 p(B|A) = \prod_{i}^{N} p(B_i|A)

と書ける*1.ただし, R_iは対戦相手のレーティングである.ここでは,http://kishi.a.la9.jp/2017R/1307.htmlの対戦当時の相手のレーティングを用いることにした.

これらを用いて計算した藤井四段のレーティングの分布 p(A|B)は次のようになる*2
f:id:kawa0616:20170816063350p:plain

また,累積分布関数は,次のようになる.
f:id:kawa0616:20170816063336p:plain

これらの結果から,藤井四段のレーティングは2030あたりである確率が最も高いことになる.
また,8/16現在レーティングが1900を超える棋士は存在しないが*3,藤井四段のレーティングが1900を超えている確率は93%ほどもあることがわかる.

結論としては,藤井四段めちゃくちゃ強い!
ただし,これは,イロレーティングの勝率の仮定が正しいときにしか当てはまらない推論であって,
この推論を他の棋士に適用したとき,実際のレーティングと整合する結果が出るのかは別途確認する必要があることに注意しておく*4

ちなみに,初戦から対局を重ねたときの(事後)確率密度関数の変化は次のようになる.
連勝中はレーティングが最大どの程度あるのか,見積もれないというのは直感に合っていて面白い.
f:id:kawa0616:20170816063624g:plain

*1:勝率の計算はElo rating system - Wikipediaを用いた.

*2:分母の積分は数値積分した

*3:棋士ランキング

*4:同様の推論を羽生三冠の直近39局でやってみたところ,レーティングが1847である確率が最も高くなった.8/16現在の羽生三冠のレーティングは1856なので,そこそこ妥当ではありそう.

各行に式番号を振った箇条書き

TeXで地味にハマったので,メモっておく.

empheqパッケージを使って実現した.

\begin{empheq}[left=\Sigma:\empheqlbrace\:]{align}
\bm x_{k+1} &=\bm A \bm x_k\label{eq:Ax}\\
y_k &= \bm c^\top \bm x_k\label{eq:Cx}
\end{empheq}

f:id:kawa0616:20170209174739p:plain

以下のサイトを参考にした.
muscle-keisuke.hatenablog.com

TeXstudio+latexmkの使い方

普段TeXstudioで原稿を書いているのですが,
日本語原稿と英語原稿で設定を切り替えたりするのが面倒なので,
latexmkを用いてディレクトリごとに設定を保持することにしました.
なお,以下の説明はMac+MacTex環境でのものです.

TeXstudioの設定

ビルド

f:id:kawa0616:20161111164536p:plain

コマンド

latexmkの項に
latexmk %.tex

latexmkrcの設置

コンパイルする文書と同じディレクトリに以下のような内容で".latexmkrc"というファイルを作ります.

日本語原稿の場合

英語原稿の場合

これで,コンパイルの設定をディレクトリごとに管理できるようになります.

ubuntu+anacondaにOpenCVをインストール(with ffmpeg)

anacondaにOpenCVをインストールしたメモ

適当にぐぐると、

conda install -c https://conda.binstar.org/menpo opencv

で簡単インストール、というのが出てくるけれど、自分の環境では動画が読み込めなかった。
どうやら、ffmpegを含めずにコンパイルしているらしい(適当)。

そこで、自前でビルドすることに。

ubuntu14.04だと、ffmpegがデフォルトではapt-get installできないようだったので、リポジトリを追加。
(これが最近までメンテされてるのかは未調査)

sudo add-apt-repository ppa:mc3man/trusty-media
sudo apt update

つぎに、
OpenCV - Community Help Wiki
にしたがって、インストールを実行。

この状態だと
/usr/local/lib/python2.7/site-packages/
以下にcv2がインストールされているため、anacondaではimportできない。

そこで、使う仮想環境のsite-packages以下にシンボリックリンクを作成する。
(参考:Install OpenCV 3 and Python 2.7+ on Ubuntu - PyimageSearch

conda create -n env27_opencv anaconda python=2.7
cd anaconda3/envs/env27_opencv/lib/python2.7/site-packages/
ln -s /usr/local/lib/python2.7/site-packages/cv2.so cv2.so

これで、anacondaの仮想環境下でcv2がimportできる。
さらに、動画の読み込みも無事成功。