弘前で働くデータサイエンティスト（？）のブログ

データ分析の勉強のための備忘録です(ノ)`ω´(ヾ)

書籍「つくりながら学ぶ! Pythonによる因果分析」の紹介

6月30日に小川雄太郎さんの新しい著作「つくりながら学ぶ! Pythonによる因果分析」が出ていたことに気づき，早速読んでみたところ，とても素晴らしい本だったので紹介記事を書くことにしました．

www.amazon.co.jp

ご存知の方も多いと思いますが，小川さんは以下の本の著者です (その他にも何冊か執筆されています)． これらの本もとても評判が良く，私もとてもお世話になりました．

今回紹介する新しい本では，因果推論や因果探索が扱われています．因果推論とは，変数間の因果関係を推定するための統計学的な方法論です．また，因果探索とは，データから因果関係の向きを推定する方法論です．ご存知の方からすれば当たり前のことですが，データから因果関係を調べたい場合，単純に回帰分析をするようなやり方では，ほとんど全ての場合に不適切です．因果推論の理論をちゃんと勉強し，適切な手法を用いる必要があります．

因果推論は一般に敷居が高いと思います．初学者にはとても分かりにくく，面倒に感じると思います．しかし，そうした方も小川さんの本なら大丈夫です．今まで因果推論をちゃんと勉強しなきゃと思いつつ逃げていた多くの方に，おすすめしたい一冊です．この本を特におすすめしたい理由は以下の4点です．

  1. 事前知識ほぼゼロの初学者でも読めること．因果推論に関する本は数学的に高度な内容を含むものが多く，全く事前知識のない方におすすめできる本を私は知りません (因果探索に関する本なら，清水昌平先生の「統計的因果探索」がとても分かりやすいですが)．一方，この本はまず，相関と因果の違いとは？という点からスタートし，初学者にも分かりやすく噛み砕いて説明してくれています．
  2. 因果推論から因果探索まで俯瞰的に学べること．伝統的な因果推論から，最近幅広い分野で応用されている因果探索まで，一冊で俯瞰的に解説した書籍は今までなかったと思います．
  3. ソースコードを自分で動かして体験できること．しかも，ただのパッケージユーザーではなく，中身が理解できるように工夫されています．
  4. 最先端の論文までカバーしていること．グラフニューラルネット，深層強化学習，GANなどを用いた因果探索手法が紹介されています．

と，ここまで書いて気づきましたが，著者の小川さん本人がQiitaに紹介記事を書かれていました．詳細はこちらをご覧下さい．

qiita.com

多次元Ornstein-Uhlenbeck過程の時間相関関数の導出メモ

応用数学

はじめに

前回の記事「多次元Ornstein-Uhlenbeck過程の確率密度関数の導出メモ」に引き続き，この記事では状態変数 \boldsymbol{x}(t)の時間相関関数を導出します．

wnk.hatenablog.com

状態変数の時間相関関数

系の定義については，前回の記事を参照して下さい．状態変数 \boldsymbol{x}(t)の時間相関関数 C(s): \mathbb{R} \to \mathbb{R}^{n \times n}を以下のように定義します．
\begin{align}
C(s) = E[\boldsymbol{x}(t+s)\boldsymbol{x}(t)^{\top}]
\end{align}
但し，ここで期待値 E[\cdot]は，定常状態 ( t \to \inftyとした場合) における期待値を意味します．このとき，この時間相関関数 C(s)は以下のように書くことができます．
\begin{align}
C(s) &= - e^{sA} {\rm vec}^{-1}( (A \oplus A)^{-1} {\rm vec}(W) ) \quad & (s \ge 0) \tag{1} \\
C(s) &= - e^{-sA^{\top}} {\rm vec}^{-1}( (A \oplus A)^{-1} {\rm vec}(W) ) \quad & (s < 0) \tag{2}
\end{align}

式(1), (2)の導出

以下のようにモーメント方程式を導くことができます．
\begin{align}
\frac{d}{ds} C(s) = \frac{d}{ds} E[\boldsymbol{x}(t+s)\boldsymbol{x}(t)^{\top}] = AC(s)
\end{align}
したがって，初期条件 C(0) = \lim_{t \to \infty} V(t) = - {\rm vec}^{-1}( (A \oplus A)^{-1} {\rm vec}(W) )の下でこの式を時間順方向に解くことにより，式(1)を得ることができます．また， C(s)の定義より C(s) = C(-s)^{\top}なので，式(2)を導くことができます．

多次元Ornstein-Uhlenbeck過程の確率密度関数の導出メモ

応用数学

はじめに

この記事では，多次元Ornstein-Uhlenbeck過程と呼ばれる確率過程について，ある初期値からこの確率過程に従って時間発展する状態変数の確率密度関数を導出します．

多次元Ornstein-Uhlenbeck過程

Ornstein-Uhlenbeck (OU) 過程は，ブラウン運動する粒子の振る舞いを記述するために提案された確率過程であり，線形な確率微分方程式によって記述されます*1．多次元の状態変数を持つOU過程を，ここでは多次元OU過程と呼び，以下の (伊藤型) 確率微分方程式で定義します．
\dot{\boldsymbol{x}}(t) = A \boldsymbol{x}(t) + \boldsymbol{w}(t) \tag{1}
但し， \boldsymbol{x}(t) \in \mathbb{R}^nはこの系の状態を表す状態変数であり，点は時間 tによる微分を表します． A \in \mathbb{R}^{n \times n}は系のダイナミクスを特徴付ける行列であり，異なる n個の固有値を持ち，全ての固有値の実部は負であるとします． \boldsymbol{w}(t) \in \mathbb{R}^nは
\begin{align}
E[\boldsymbol{w}(t)] = 0, \ E[\boldsymbol{w}(t)\boldsymbol{w}(s)^{\top}] = W\delta(t-s)
\end{align}
を満たす白色ガウスノイズとします ( \delta(\cdot)はDiracのデルタ関数です)．

状態変数の確率密度関数

状態変数 \boldsymbol{x}(t)がある初期値 \boldsymbol{x}(0)から式(1)に従って時間発展するとします．このとき，時刻 t \ge 0の状態変数 \boldsymbol{x}(t)は，ガウス分布の再生性より多次元ガウス分布に従います．この多次元ガウス分布の平均 \boldsymbol{m}(t)，共分散行列 V(t)を以下のように定義します．
\begin{align}
\boldsymbol{m}(t) &= E[\boldsymbol{x}(t)] \\
V(t) &= E[(\boldsymbol{x}(t) - \boldsymbol{m}(t)) (\boldsymbol{x}(t) - \boldsymbol{m}(t))^{\top}]
\end{align}
このとき， \boldsymbol{m}(t)と V(t)の時間発展は以下のように記述されます．
\begin{align}
\boldsymbol{m}(t) &= e^{tA} \boldsymbol{x}(0) \tag{2} \\
V(t) &= {\rm vec}^{-1}\left( \int_0^t e^{s(A \oplus A)} {\rm vec}(W) ds \right) \tag{3}
\end{align}
但し， {\rm vec}(\cdot):\ \mathbb{R}^{n \times n} \to \mathbb{R}^{n^2}はある行列 M=(M_{i,j}) \in \mathbb{R}^{n \times n}について
\begin{align}
{\rm vec}(M) = [M_{1,1}, \ldots, M_{n,1}, M_{1,2}, \ldots, M_{n,2}, \ldots, M_{1,n}, \ldots, M_{n,n}]^{\top}
\end{align}
と定義される関数であり， \oplusはKronecker和*2を表します．また特に， t \to \inftyのとき，
\begin{align}
\boldsymbol{m}(t) &\to \boldsymbol{0} \\
V(t) &\to -{\rm vec}^{-1}\left( (A \oplus A)^{-1} {\rm vec}(W) \right)
\end{align}
となります．

式(2), (3)の導出

以下のように， \boldsymbol{m}(t)と V(t)の時間発展を記述する方程式を得ることができます．なお，式(5)の導出には伊藤の補題*3を利用します．
\begin{align}
\dot{\boldsymbol{m}}(t) &= \frac{d}{dt} E[\boldsymbol{x}(t)] = A \boldsymbol{m}(t) \tag{4} \\
\dot{V}(t) &= \frac{d}{dt} E[(\boldsymbol{x}(t) - \boldsymbol{m}(t)) (\boldsymbol{x}(t) - \boldsymbol{m}(t))^{\top}] = AV(t) + V(t)A^{\top} + W \tag{5}
\end{align}
式(4)を初期条件 \boldsymbol{m}(0) = \boldsymbol{x}(0)の下で解くことにより，式(2)を導くことができます．また，式(5)はKronecker和を用いることにより，以下のように書き直すことができます．
\begin{align}
\frac{d}{dt} {\rm vec}(V(t)) = (A \oplus A) {\rm vec}(V(t)) + {\rm vec}(W)
\end{align}
この式を初期条件 V(t) = Oの下で解くことにより，式(3)を導くことができます．

*1:Stochastic dynamical systems - Scholarpedia

*2:Kronecker product - Wikipedia

*3:Itô's lemma - Wikipedia

青森市周辺の観光スポットを紹介

雑記

私の結婚式で青森市を訪れる友人たち向けに、青森市周辺の観光地をまとめました。 青森は私たちの大好きな場所なので、ぜひ楽しんで行ってもらえればと思います。 青森といえば、やはり魚介類と温泉です。

青森市街地

戦災により歴史的な建造物はほとんど残っていません。 なので、あまり観光向きではありませんが、青森の新鮮な魚介類を楽しむことができます。

弘前市街地 (青森からJR奥羽線で40分)

弘前は、お城とその城下町や、明治期の洋館が残る美しい街です。 漫画・アニメ「ふらいんぐうぃっち」の聖地でもあります。 http://www.hirosaki-kanko.or.jp/web/edit.html?id=fc_flyingwitch

弘前・岩木山 (青森から自動車で1.5時間)

霊峰・岩木山は、古くから津軽の人々の信仰の対象でした。 美しい景色、極上の釜飯に、最高の温泉をお楽しみ下さい。 (※レンタカーの利用を推奨します)

八甲田山・奥入瀬渓流・十和田湖 (青森からバスで1〜1.5時間)

青森観光のメインストリートです。 青森市街地から南へ山道を進むと、順に八甲田山、奥入瀬渓流、十和田湖があります。 (※バス移動も可能ですがレンタカーがあると便利です)

バス情報: http://www.jrbustohoku.co.jp/route/detail/?RID=1

以上、青森市街地から1時間程度で行くことのできる範囲をざっと紹介しました。より詳しいことについても、もしよければご相談に乗ります。是非、青森を楽しんで行って下さい！

追記(2018/6/11): @yamas4kiさんや@takatoh1さんに「浅虫温泉が入っていない！」というお叱り(?)をいただいたので追記します。

夏泊半島 (青森から青い森鉄道で30分)