何も解らない。。

データサイエンス関連

カルマンフィルター

【カルマンフィルター】

観測方程式: \begin{equation} \begin{split} Y_{t} = \beta_{t}X_{t} + e_{t} \end{split} \end{equation}

状態方程式: \begin{equation} \begin{split} \beta_{t} = T\beta_{t-1} + \epsilon_{t} \end{split} \end{equation}

この時、カルマンゲインを \begin{equation} \begin{split} K_{t} := \frac{X_{t}Σ_{t|t-1}}{X_{t}^2Σ_{t|t-1}+σ_{e}^2} \end{split} \end{equation} と表すと、状態変数の期待値、分散のフィルタリングは以下、

\begin{equation} \begin{split} \hat{\beta_{t|t}} = \hat{\beta_{t|t-1}} + K_{t}(y_{t} - \hat{Y_{t|t-1}}) \cr \hat{Σ_{t|t}} = (1 - X_{t}K_{t}) \hat{Σ_{t|t-1}} \end{split} \end{equation}

で表せる。

適当な紹介にも程があります。順を追ってみていきましょう。

今回は以下の本を参考にまとめます。

経済・ファイナンスのための カルマンフィルター入門 (統計ライブラリー)

経済・ファイナンスのための カルマンフィルター入門 (統計ライブラリー)

Amazon

はじめに

本記事の構成ですが、まず

  1. 状態空間モデルについて簡単に説明
  2. フィルタリングの直感的理解
  3. フィルタリングの公式の導出

という流れになっています。

尚、上記のフィルタリングの式の導出では、簡単のため、幾つかの仮定*1をおいたものになっております。あしからず。

そもそも状態空間モデルってなんだっけ

観測した結果についての関係式”観測方程式”と、観測できない状態に関する関係式”状態方程式”からなる時系列分析モデル、です。

…と言われてもなかなかピンと来ませんよね。。2つの式を再掲します。

観測方程式： Y_{t} = \beta_{t}X_{t} + e_{t}

状態方程式： \beta_{t} = T\beta_{t-1} + \epsilon_{t}

具体的に考えてみます。

例えば、消費関数*2を例に考えてみます。 Y_{t}が消費、 X_{t}が所得、 \beta_{t}を限界消費性向とします。よく見る式だと上の観測方程式オンリーで \beta_{t}が定数 \betaになったものではないでしょうか。

Y_{t} = \beta X_{t} + e_{t}

つまり状態方程式のないノーマルな消費関数の式*3では限界消費性向が時間で不変であると仮定してモデリングをしていることになります。

これを状態方程式を考えることで、観測できない状態（限界消費性向）を時間で変化するものとして扱うことができます。このように状態空間モデルを使うことで、より柔軟なモデリングが可能になるというわけです。

フィルタリングとは（直感的な理解）

カルマンフィルターの目的は、確率変動する状態変数 \beta_{t}の平均と分散を推定することにあります。

本記事ではそのような推定の一つである”フィルタリング”をメインに解説しますが、他にも”一期先予測”と”スムージング”という推定があります。これらの直感的な関係性を参考書籍に出てくる例で確認します。本の中ではこれらの推定を”ミステリー小説の犯人予想”に喩えて説明しています。

まず一期先予測ですが、これは分かりやすい推定で、例えば”99ページまで本を読んだ状態で100ページ目で誰が犯人か予測する”というものです。 次にスムージングは、本を全て読み終えた状態で100ページ目*4での犯人予測がなぜ間違いだったか（正解だったか）を振り返るような作業にあたります。最後に本記事のメインのフィルタリングは、100ページ目で得た情報をもとに一期先予測の修正を行うような作業、と考えればなんとなく理解ができると思います。

ベイジアンアプローチによるフィルタリング公式の導出

さて、登場人物の説明が概ね完了したので、最後に実際にどう推定（フィルタリング）するのか、というところを詰めていきましょう。

参考書籍では、正規性を仮定し導出、正規性を仮定せず平均分散から導出、ベイズの定理から導出、の3タイプの導出が掲載されています。*5 本記事では、この中から”ベイズの定理から導出”（ベイジアンアプローチ）を紹介することにします。

前提条件の確認

まず、前提条件について確認しましょう。観測、状態方程式の誤差項について以下を仮定します。

\begin{equation} \begin{split} e_{t} \propto N(0, σ_{e}^2), \epsilon_{t} \propto N(0, σ_{\epsilon}^2) \cr E[e_{t}\epsilon_{s}] = 0 \qquad \forall s, t \cr E[e_{t}e_{s}] = 0, E[\epsilon_{t}\epsilon_{s}] = 0 \quad \forall s \neq t \end{split} \end{equation}

誤差項は正規分布に従い、系列相関なし、観測、状態方程式のそれぞれの誤差同士も相関を持たないということですね。 この前提のもと、時点 t-1で時点 tの予測を行って、時点 tで得られた観測値Y_{t}をもとに予測値の更新を行う…という作業を逐次的に行っていくことになります。

ベイズの式とカルマンフィルターの関係

次に今回の記事の主役であるベイズの式とカルマンフィルターの関係性を考えると以下のようになります。

Pr(\beta_{t}|Y_{t}=y_{t})\propto Pr(Y_{t}|\beta_{t})Pr(\beta_{t})

カルマンフィルタの文脈で言うと、事前に予測していた状態変数（事前分布）を新たに観測した情報（尤度）をもとに現時点での状態変数の予測（事後分布）に更新（フィルタリング）している式だと読み取れます。

それでは事前分布と尤度を具体的に見ていきます。

事前分布

事前分布 Pr(\beta_{t})は、 Y_{t} = yが観測される前の時点（ t-1）での \beta_{t}の分布です。

期待値と分散は時点 t-1で活用できる情報を使って推定される値、 \hat{\beta_{t|t-1}}, \hat{Σ_{t|t-1}}で、仮定よりこの分布は正規分布に従うため、 \hat{\beta_{t}}|\Omega_{t-1} \sim N(\hat{\beta_{t|t-1}}, \hat{Σ_{t|t-1}})となります。

事前分布： Pr(\beta_{t}) = \frac{1}{\sqrt{2\pi \hat{Σ_{t|t-1}}}} exp(-\frac{1}{2} \frac{ {(\beta_{t} - \hat{\beta_{t|t-1}})}^2 }{\hat{Σ_{t|t-1}} })

尤度

尤度 Pr(Y_{t}|\beta_{t})は、 \beta_{t}がわかった状態で観測変数 Y_{t}を考えます。

観測方程式： Y_{t} = \beta_{t}X_{t} + e_{t}

であることを思い出すと、尤度は Y_{t}|\beta_{t} \sim N(\beta_{t}X_{t}, σ_{e}^2)となります。

尤度： Pr(Y_{t}|\beta_{t}) = \frac{1}{\sqrt{2\pi σ_{e}^2}} exp(-\frac{1}{2} \frac{ {(Y_{t} - \beta_{t}X_{t})}^2 }{ σ_{e}^2 })

事後分布

最後に事後分布ですが、事後分布は2通りの方法で表現出来ます。

まずは今まで見てきた事前分布、尤度を使って以下のように表現出来ます。

\begin{equation} \begin{split} Pr(\beta_{t}|Y_{t}) \propto \frac{1}{\sqrt{2\pi σ_{e}^2}} exp(-\frac{1}{2} \frac{ {(Y_{t} - \beta_{t}X_{t})}^2 }{ σ_{e}^2 }) \times \frac{1}{\sqrt{2\pi \hat{Σ_{t|t-1}}}} exp(-\frac{1}{2} \frac{ {(\beta_{t} - \hat{\beta_{t|t-1}})}^2 }{\hat{Σ_{t|t-1}} }) \cr \propto exp(-\frac{1}{2}(\frac{X_{t}^2}{σ_{e}^2} + \frac{1}{Σ_{t|t-1}})\beta_{t}^2 + (\frac{X_{t}y_{t}}{σ_{e}^2} + \frac{\beta_{t|t-1}}{Σ_{t|t-1}})\beta_{t}) \end{split} \end{equation}

細かい計算はキャッツアイしますが、指数部分を \betaに関して整理していけばこのようになります。

もう一つの方法は、

条件付き期待値： \hat{\beta_{t|t}} \equiv E[\beta_{t}|Y_{t} = y_{t}]

分散： \hat{Σ_{t|t}} \equiv Var[\beta_{t}|Y_{t} = y_{t}]

を使って表現します。

\begin{equation} \begin{split} Pr(\beta_{t}|Y_{t}) = \frac{1}{\sqrt{2\pi \hat{Σ_{t|t}}}} exp(-\frac{1}{2} \frac{ {(\beta_{t} - \hat{\beta_{t|t}})}^2 }{ \hat{Σ_{t|t}^2} }) \cr \propto exp(-\frac{1}{2}\frac{\beta_{t}^2 -2\beta_{t}\hat{\beta_{t|t}}}{\hat{Σ_{t|t}}}) \cr = exp(-\frac{1}{2} \frac{1}{\hat{Σ_{t|t}}} \beta_{t}^2 + \frac{\hat{\beta_{t|t}}}{\hat{Σ_{t|t}}} \beta_{t}) \end{split} \end{equation}

こちらも粛々と計算するのみですが、ポイントとしては確率変数 \beta_{t}に注目して式変形していくところでしょうか。推定値である、 \hat{\beta_{t|t}} 、 \hat{Σ_{t|t}}は無視します。

これで準備が整いました。あとは上記の2式を比較することで、フィルタリングの公式を導出しましょう。

フィルタリングの公式

上記の事後分布の2通りの表現について、 \beta_{t}と \beta_{t}^2に注目して比較することで、以下の式を得ます。

\begin{equation} \begin{split} \frac{1}{\hat{Σ_{t|t}}} = \frac{X_{t}^2}{σ_{e}^2} + \frac{1}{\hat{Σ_{t|t-1}}} \cr \frac{\hat{\beta_{t|t}}}{\hat{Σ_{t|t}}} = \frac{X_{t}y_{t}}{σ_{e}^2} + \frac{\hat{\beta_{t|t-1}}}{\hat{Σ_{t|t-1}}} \end{split} \end{equation}

これを式変形することで冒頭のフィルタリングの式が導出できます。まず、分散のフィルタリングの式から導出しましょう。

上記の最初の式を変形することで以下を得ます。

\begin{equation} \begin{split} \hat{Σ_{t|t}} = \frac{σ_{e}^2}{X_{t}^2 \hat{Σ_{t|t-1}} + σ_{e}^2} \hat{Σ_{t|t-1}} \cr = (1 - X_{t}\frac{X_{t}\hat{Σ_{t|t-1}}}{X_{t}^2 \hat{Σ_{t|t-1}} + σ_{e}^2}) \hat{Σ_{t|t-1}} \cr = (1 - X_{t}K_{t}) \hat{Σ_{t|t-1}} \end{split} \end{equation}

1行目から2行目はややテクニカルですが、1行目の分子に X_{t}^2 \hat{Σ_{t|t-1}} - X_{t}^2 \hat{Σ_{t|t-1}} = 0を足してみると分かりやすいかもしれません。

最後に期待値のフィルタリングの式を導出します。こちらは少し厄介ですが、頑張って計算します。

\begin{equation} \begin{split} \hat{\beta_{t|t}} = \frac{X_{t} \hat{Σ_{t|t}}}{σ_{e}^2} y_{t} + \frac{\hat{Σ_{t|t}}}{\hat{Σ_{t|t-1}}} \hat{\beta_{t|t-1}} \cr = \hat{Σ_{t|t}} (\frac{X_{t} y_{t} \hat{Σ_{t|t-1}} + \hat{\beta_{t|t-1}} σ_{e}^2}{σ_{e}^2 \hat{Σ_{t|t-1}}}) \cr = \frac{σ_{e}^2}{X_{t}^2 \hat{Σ_{t|t-1}} + σ_{e}^2} \hat{Σ_{t|t-1}} (\frac{X_{t} y_{t} \hat{Σ_{t|t-1}} + \hat{\beta_{t|t-1}} σ_{e}^2}{σ_{e}^2 \hat{Σ_{t|t-1}}}) \cr = \frac{\hat{\beta_{t|t-1}} σ_{e}^2}{X_{t}^2 \hat{Σ_{t|t-1}} + σ_{e}^2} + \frac{X_{t} \hat{Σ_{t|t-1}}}{X_{t}^2 \hat{Σ_{t|t-1}} + σ_{e}^2} y_{t} \cr = \frac{\hat{\beta_{t|t-1}} σ_{e}^2}{X_{t}^2 \hat{Σ_{t|t-1}} + σ_{e}^2} + K_{t}y_{t} + (K_{t}\hat{\beta_{t|t-1}}X_{t} - K_{t}\hat{\beta_{t|t-1}}X_{t}) \cr = \frac{σ_{e}^2}{X_{t}^2 \hat{Σ_{t|t-1}} + σ_{e}^2} \hat{\beta_{t|t-1}} + \frac{X_{t} \hat{Σ_{t|t-1}}}{X_{t}^2 \hat{Σ_{t|t-1}} + σ_{e}^2} X_{t} \hat{\beta_{t|t-1}} + K_{t}(y_{t} - \hat{\beta_{t|t-1}} X_{t}) \cr = \hat{\beta_{t|t-1}} + K_{t}(y_{t} - \hat{Y_{t|t-1}}) \end{split} \end{equation}

こちらの式変形でも K_{t}\hat{\beta_{t|t-1}}X_{t} - K_{t}\hat{\beta_{t|t-1}}X_{t} = 0をうまく使って式変形しています。

以上で目的のフィルタリング式が導出できました。お疲れ様でした。

記念に式を再掲しておきましょう。

\begin{equation} \begin{split} \hat{\beta_{t|t}} = \hat{\beta_{t|t-1}} + K_{t}(y_{t} - \hat{Y_{t|t-1}}) \cr \hat{Σ_{t|t}} = (1 - X_{t}K_{t}) \hat{Σ_{t|t-1}} \end{split} \end{equation}

なんだか愛着が湧いてきましたね。カルマンフィルターは個人的には小難しくてなんだかとっつきにくい印象だったのですが、この参考書籍のおかげで少し仲良くなれた気がします。 具体例なども豊富に取り扱っていて、読むためのバーもそんなに高くないので、とてもおすすめできる逸品です。

それでは今回はこの辺で、ごきげんよう。

*1:具体的には、状態変数、観測変数がそれぞれ1つで、正規分布に従うことを仮定しています。正規性を仮定しないベイジアンアプローチについては、谷崎（1993）の第１章、北川（2005）の第14章を参照、とのことです。

*2:上で紹介した本でも事例として出てきます。

*3:別に限界消費性向が定数なものがデファクトというわけではないでしょうけど。。

*4:特に100ページ目でなくても良いです。

*5:とても贅沢な内容になっています。

条件付き期待値

【条件付き期待値】

E[Y|X] = E[Y|X = x] = \int yf_{Y|X=x}(y) dy

今回はとてもシンプルですね。（というか定義を載せただけ） ただこの概念、初学だと結構混乱すると思うんですよね。。

特に条件付き期待値の期待値を考えたり、条件付き分散を考えたり、、「なんとなく」で理解していると結構躓く概念だと思うので、式の意味を考えたり、条件付き期待値にまつわる有名な公式を証明したりしながら、苦手意識のある方の症状緩和を目標にやっていこうと思います。

まず、条件付き”確率（密度関数）”を改めて見てみましょう。

f_{Y|X}(y) = \frac{f(x, y)}{f_{X}(x)}

ここで大事なのが、 f_{Y|X}(y)が yの関数だということです。

というのもこれから見ていく条件付き確率の公式もそうですが、 Y|Xの Yも Xもどちらも動かして考えるので、今動かして考えているのはどちらなのか、しっかり意識して式を追わないと迷子になりやすい*1概念だと思っています。

上の条件付き確率では、 yの関数であることから動かすのは Yの方、 X = xで固定した（条件づけられた） Yの分布を考えているということになります。

条件付き期待値の話に戻ります。条件付き期待値 E[Y|X]は先ほど確認したとおり、 Yが主役（？）です。そう考えると、 \int yf_{Y|X=x}(y) dyで定義されるのも違和感はないんじゃないでしょうか。

では、準備運動はここまでとして、次のやや込み入った話題に移りましょう。

【期待値の繰り返しの公式（全確率の公式）】

E[E[Y|X]] = E[Y]

これ、混乱しませんか。。でも大丈夫です。先ほどの話をもとにどの変数を動かしているのか明示的に式に書くことで、頭を整理しましょう。

E_{X}[E_{Y}[Y|X]] = E_{Y}[Y]

これでハッキリしました。一応式を追っておくと、中身の条件付き期待値が Xで固定して Yの期待値を取ってます。で、その結果に対して今度は Xを動かして期待値を取ってます。ここいらで証明をしておきましょう。

【証明】

\begin{equation} \begin{split} E_{X}[E_{Y}[Y|X]] = \int E_{Y}[Y|X]f_{X}(x)dx \cr = \int (\int y\frac{f(x, y)}{f_{X}(x)}dy) f_{X}(x) dx \cr = \int \int yf(x, y) dx dy \cr = E_{Y}[Y] \end{split} \end{equation}

結構シンプルに証明できました。ではこの流れで条件付き分散を使った公式も見ていきましょう。 まず条件付き分散の定義です。

【条件付き分散】

Var_{Y}[Y|X] = E_{Y}[ (Y - E_{Y}[Y|X ] )^2|X]

今までの要領で式を眺めれば良いので、これは特に説明は不要でしょう。では、公式を見ていきます。

Var_{Y}[Y] = E_{X}[Var_{Y}[Y|X]] + Var_{X}[E_{Y}[Y|X]]

Yの分散は、Yの条件付き分散の期待値 + Yの条件付き期待値の分散と分解できる、 ということになります。本当でしょうか、証明を追ってみましょう。

【証明】

\begin{equation} \begin{split} Var_{Y}[Y] = E_{Y}[(Y-E_{Y}[Y])^2] \cr = E_{X} [ E_{Y}[(Y-E_{Y}[Y|X]+E_{Y}[Y|X]-E_{Y}[Y])^2 ] | X ] \cr = E_{X} [ E_{Y}[(Y-E_{Y}[Y|X])^2 | X ]] + E_{X} [ E_{Y}[(E_{Y}[Y|X]-E_{Y}[Y])^2 ] | X ]] \cr + 2E_{X} [ E_{Y} [(Y-E_{Y}[Y|X])(E_{Y}[Y|X]-E_{Y}[Y]) | X ]] \cr = E_{X} [ Var_{Y}[Y|X] ] + Var_{X} [ E_{Y}[ Y|X ]] \end{split} \end{equation}

2行目でちゃっかり期待値の繰り返しの公式を使っているのがポイントですね。

最後の式変形は少し補足が必要かもしれません。

まず、 2E_{X}の項が消えているのは、

\begin{equation} \begin{split} 2E_{X} [ E_{Y} [(Y-E_{Y}[Y|X])(E_{Y}[Y|X]-E_{Y}[Y]) | X ]] \cr = 2E_{X} [(E_{Y}[Y|X]-E_{Y}[Y]) E_{Y} [Y-E_{Y}[Y|X] | X ]] \cr = 0 \cr \cr \because E_{Y} [Y-E_{Y}[Y|X] | X ] = 0 \end{split} \end{equation}

ここの処理は今回の記事の集大成のような感じですね。 Eの影響する範囲はどこか、（この Eは、 E_{x} , E_{y}のどちらなのか）意識して計算すればややこしいですが追えますね。

残った部分は定義に従って変形しているだけですが、2項目は、

E_{X}[E_{Y}[Y|X]] = E_{Y}[Y]

となることに注意して分散の形に変形します。

はい、これで当初目標としていた分散の公式まで証明できました。かなり冗長に書いたので逆に見にくいところもありますが、実際に式を追うときは、条件付き期待値の値などを適宜記号などで置き換えるともう少し簡潔にまとまると思います。

この公式、今回はちゃんと証明しましたが、やや煩雑なので一回導出したら覚えちゃってもいいかもしれません。*2

Yの分散は、 Y|Xの EVと VEの和と等しい。意外と覚えやすいと思います。もちろん Eと Vの添字（ Xか Yか）はどっちがどっちかは、もう迷わないようになりましたよね？

では今回はこの辺で、そりでわ。

*1:と少なくとも私は感じました。。

*2:暇なら年に一回くらい証明しましょう。

ネイマン・ピアソンの補題

【ネイマン・ピアソンの補題】

f(x, \theta_{i}), i = 0,1を帰無仮説及び対立仮説のもとでの密度関数（あるいは確率関数）とする。与えられた c \ge 0と r \ (0 \le r \le 1)に対してサイズが \alphaである次のような検定関数を考える。

\begin{align} \delta_{c,r}(x) = \begin{cases} 1 & (\frac{f(x, \theta_1)}{f(x, \theta_0)} > c)\cr r & (\frac{f(x, \theta_1)}{f(x, \theta_0)} = c)\cr 0 & (\frac{f(x, \theta_1)}{f(x, \theta_0)} < c) \end{cases} \end{align}

このとき、有意水準 \alphaの検定のなかで、 \delta_{c,r}が最強力検定となる。

はい、いきなりそんなこと言われても困ります。ということで、毎度お馴染み弊家本棚の積読筆頭候補の『現代数理統計学』を参考に順を追ってみていきましょう。

新装改訂版 現代数理統計学

新装改訂版 現代数理統計学

Amazon

まず、最強力検定とは何ぞやという話ですが、その前に少し準備をします。

検定関数：決定関数（ \delta : X \to D）

帰無仮説を棄却するなら1、受容するなら0を取る関数。

リスク関数： 損失関数の期待値（ R(\theta, \delta) = E_{\theta}[L(\theta, \delta(X))]）

損失関数は、正しく判断ができなかった時に1をとる0-1損失関数とします。 期待値を取っているのは、 \thetaが未知であり、損失関数の値自体も確率的に変動するので、単純に損失を最小化しよう、とはできないわけです。そのため平均的な損失を考えるか、ということでこのリスク関数というものが定義されています。

検出力関数： \beta_{\delta}(\theta) = E_{\theta} [\delta(X)] = P_{\theta}(\delta(X) = 1)

\beta_{\delta}(\theta)を検出力関数といいます。 尚、先に見たリスク関数は、検出力関数を使って以下のように書けます。

\begin{align} R(\theta, \delta) = \begin{cases} \beta_{\delta}(\theta) & (\theta \in \Theta_{0})\cr 1- \beta_{\delta}(\theta) & (\theta \in \Theta_{1}) \end{cases} \end{align}

はい、話を戻して最強力検定について確認していきます。

\deltaを任意の有意水準 \alphaの検定、すなわち \beta_{\delta}(\theta) \le \alpha, \ \forall \theta \in \Theta_{0}とする。

この時、 \delta^*が有意水準 \alphaの一様最強力検定（UMP）であるとは、

\beta_{\delta^*}(\theta) \ge \beta_{\delta}(\theta), \ \forall \theta \in \Theta_{1} が成立する事である。

ここで、”一様”というのは全ての対立仮説（ \forall \theta \in \Theta）について検出力を最大化する、という意味であり、単純仮説（ \Theta_{1} = \theta_{1}）の場合は単に最強力検定と呼ばれるようです。

一旦ここまでの話を整理して、もう一度ネイマン・ピアソンの補題を見てみましょう。

サイズ（第1種の過誤の上限）が \alphaであるような検定 \delta_{c, r}を考えます。その中で最も検出力の高い検定を最強力検定と言います。ところで、冒頭で挙げたような検定関数を構成すると、そいつが噂の最強力検定です。

なんとか言ってることは分かりました。あとはこれを証明しましょう。 連続分布の場合の証明になります。（離散の場合も同様の議論）

【証明】

\deltaを有意水準 \alphaの検定関数とした時、

\int (\delta_{c,r}(x) - \delta(x))(f(x, \theta_{1}) - cf(x, \theta_{0})) dx \ge 0 \dots (*) が成立する（これは後で示す）が、これを以下のように整理することで題意が示せる。

\begin{equation} \begin{split} \beta_{\delta_{c,r}}(\theta_{1}) - \beta_{\delta}(\theta_{1}) = \int (\delta_{c,r}(x) - \delta(x))f(x, \theta_{1})dx \cr \ge c(\int (\delta_{c,r}(x) - \delta(x))f(x, \theta_{0})dx) \cr = c(E_{\theta_{0}} [\delta_{c,r}(X)] - E_{\theta_{0}} [\delta(X)] ) \cr = c(\alpha - E_{\theta_{0}} [\delta(X)] ) \cr \ge 0 \end{split} \end{equation}

(*)の証明

0 \le \delta(x) \le 1であることに注意して、

(1) \ f(x, \theta_{1}) > cf(x, \theta_{0})の時

\delta_{c,r}(x) = 1（そのように検定関数を作ったのでした）であることと、 0 \le \delta(x) \le 1、また条件の f(x, \theta_{1}) > cf(x, \theta_{0})を合わせて考えて、 (*)の被積分関数は非負になることが分かる。

(2) \ f(x, \theta_{1}) ＜ cf(x, \theta_{0})の時

これも同様の議論により (*)の被積分関数は非負になることが分かる。

(3) \ f(x, \theta_{1}) = cf(x, \theta_{0})の時

この場合、 (*)の被積分関数は0になる。

以上の結果から (*)が示せた。

はい、証明まで一通り追ってみました。*1 ネイマン・ピアソンと少しは仲良くなれそうでしょうか。

では今回はこの辺で、そりでわ。

*1:どうでもいいですけど、tex書くのってなかなか時間溶かしますね。。

十分統計量

【十分統計量】

k個の統計量 T = (T_{1}, ..., T_{k})がパラメータθに関するk次元の十分統計量であるとは、Tを与えたときの X = (X_{1}, ..., X_{n})の条件つき分布が \thetaに依存しないことである。


初見だと？ってなりませんか、私はなりました。
ということで、この機会に弊家本棚の積読筆頭候補の『現代数理統計学』を参考にまとめてみました。

新装改訂版 現代数理統計学

新装改訂版 現代数理統計学

Amazon

要は、 \thetaで特徴付けられる分布があるんだけど、十分統計量Tなるものの情報を得るとその条件つき分布は \thetaに依存しなくなる、ということのようです。

まだピンと来ないので具体例を見ていきたいと思いますが、その前に『分解定理』を紹介して見通しを良くしましょう。

【分解定理】

Xを離散確率変数または連続確率変数とし p_{\theta}をXの確率変数または密度関数とする。 T(X) = (T_{1}(X), ..., T_{k}(X))が十分統計量であるための必要十分条件は p_{\theta}(x)が 、


p_{\theta}(x) = g_{\theta}(T(x))h(x)


の形に分解できることである。
ここで h(x)は \thetaを含まないxのみの関数である。


密度関数を \thetaを含むxの統計量Tの関数と含まないxの関数に分けて、前者の統計量が \thetaの十分統計量になっている、と。

この定理の嬉しいところは、Tが十分統計量かどうかを調べるために条件付き分布（ P_{\theta}(X = x | T = t)）を求める必要がない、という点にあります。特に連続変数の場合は条件付き分布の導出に（自明でない）変数変換が必要となり、難儀するとのこと。

この定理は連続変数版の証明は測度論が必要とのことなので、離散版だけ証明を追いましょう。

【証明】

（必要性）
p_{\theta}(x) = g_{\theta}(T(x))h(x)の分解が出来ているとする。


p_{\theta}(T = t) = \sum_{T(x)=t} p_{\theta}(x) = \sum_{T(x)=t} g_{\theta}(T(x))h(x) = g_{\theta}(t) \sum_{T(x)=t}h(x)


なので、


p_{\theta}(X = x | T = t) = \frac{g_{\theta}(t)h(x)}{g_{\theta}(t)\sum_{T(x)=t} h(x)} = \frac{h(x)}{\sum_{T(x)=t} h(x)}


となり、 \thetaに依存しないことから、Tが十分統計量となる。


（十分性）
Tが十分統計量とすると、


p_{\theta}(T = t) = g_{\theta}(t), p_{\theta}(X = x | T = t) = h(x)（ \thetaに依存しない）と置いて、


p_{\theta}(X = x) = p_{\theta}(T = t) * p(X = x | T = t) = g_{\theta}(T(x))h(x)


となる。


さて、話を戻してポアソン分布（ X_{1}, ..., X_{n} \sim Po(\lambda)）について具体例を見ていきましょう。

p_{\lambda}(x) = \prod_{i = 1}^n \frac{\lambda^x_{i}}{x_{i}!} e^{-\lambda} = \lambda^{\sum_{i = 1}^n}e^{-n\lambda}(\prod_{i = 1}^n x_{i}!)^{-1}


パラメータ \lambdaの入った部分とそれ以外に分けて、 g_{\lambda} = \lambda^{\sum_{i = 1}^n}e^{-n\lambda}, h(x) = (\prod_{i = 1}^n x_{i}!)^{-1}と置けば、分離定理より \sum_{1}^n X_{i}は十分統計量だと分かりました。

こうして具体例を見ると大分クリアになりますね。

せっかくだから6章全部まとめようと思ったんですが、前章までの知識が必要だったりしてすぐ纏まんなそうなので、今回はこの辺で、そりでわ。

『RとStanではじめるベイズ統計モデリングによるデータ分析入門』

実践Ｄａｔａ　Ｓｃｉｅｎｃｅシリーズ　ＲとＳｔａｎではじめる　ベイズ統計モデリングによるデータ分析入門 (ＫＳ情報科学専門書)

実践Ｄａｔａ　Ｓｃｉｅｎｃｅシリーズ　ＲとＳｔａｎではじめる　ベイズ統計モデリングによるデータ分析入門 (ＫＳ情報科学専門書)

統計モデリング系のお勉強は、昔『みどり本』*1を読んだ後『アヒル本』*2を半分くらい読んでほったらかしにしていました。実務で必要に迫られたこともなかったので、いつかやろうと先延ばしにしていたのですが、今回良い機会なので、この本から再入門することにしました。

勉強した本の紹介方法は、現状手探りでやっていこうかなと思っていて、今回は各部毎の内容を一部抜粋でまとめつつ、最後に全体の感想を書くようなスタイルでやってみようかなと思います。*3

尚、私の読解力、知識不足により誤読している箇所があるかもしれないので、もし見つけられた方は優しく諭して頂けますと幸いです。

内容紹介

1部 【理論編】ベイズ統計モデリングの基本

1部は理論編ということで、確率や統計の基礎事項、モデリングとはなんぞや、事後分布に従う乱数を作成するためのMCMC、などをコンパクトに紹介した部です。

”5章からはちゃんと読んでね”とある通り、ここから統計モデル構築のための準備が始まります。大事なのがベイズ更新の考え方ですかね。

ベイズ更新

P(H_1|D) = P(H_1)\dfrac{P(D|H_1)}{P(D)}

P(H_1)が事前確率、つまりデータが得られる前に想定していた確率ですね。 P(H_1|D)が事後確率、つまりデータが得られた後に想定する確率となります。これらを結びつけるのが、\dfrac{P(D|H_1)}{P(D)}の部分です。分子が尤度、分母が周辺尤度となります。つまり考えうる様々な仮定の元でデータが得られる確率のうち、今回想定している仮定H_1の元でデータが得られる確率を考えることで、事前分布を事後分布に更新しているわけです。

ベイズ推論とMCMC

パラメータの事後分布が推定できれば、万事解決のように見えるんですが、まだ問題があります。それが、事後分布が複雑過ぎて扱いにくい問題です。例えばパラメータθの期待値を求めるのに

\int_{-∞}^{∞}f(\theta｜D)\cdot\theta d\theta

この計算をしないといけないですが、先程の事後分布が複雑過ぎ問題があるため容易に計算ができなかったりします。そこで、事後分布に従う乱数をいくつか生成して、その平均を取ることで代替しようということを考えます。そのために使うのがMCMCというわけですね。

乱数発生アルゴリズムについては、本書ではMH法とHMC法が紹介されています。MH法のざっくりした考え方としては、以下のような感じです。

  1. パラメータの初期値を設定する
  2. パラメータの値を少しズラす
  3. ズラす前後でカーネルを計算する
  4. カーネルの比をとって、値を更新するか決める

手順の2.~3.を適当な回数繰り返して乱数を得ます。カーネルというのは、事後分布の分子であるP(\theta)P(D|\theta)のことです。分母の周辺尤度は無くて良いのかと思われるかもしれませんが、これは上記の手続きで”比”をとっているため必要がないのですね。この点は結構大事で、そもそも周辺尤度の計算自体が難儀なため、使わないで済むととても有り難いわけです。

2部 【基礎編】RとStanによるデータ分析

2部は前半は基本的なR、可視化周り、stanの使い方の解説が中心です。後半では実際にstanを使って簡単な事例に取り組みます。そうそう、こんな感じだったなあとか言いながら写経しました。対数密度加算文形式の表記もちゃんと紹介されています。

target += normal_lpdf(データ|パラメータ)

3部 【実践編】一般化線形モデル

いよいよ実践編です。まずは一般化線形モデルについて概要がまとめられており、その後実際にStanによる分析例をいくつか追っていくような形式です。こういうのは慣れもあるので、適度な分量を写経できる感じになってて有り難いです。

最後の交互作用の章にある売り上げと製品数、店員数の関係性の分析は、交互作用を考えることで関係性が浮き彫りになる面白い例です。また、結果が因果関係ではないことを指摘されているのもとても大事なポイントだと思います。*4

線形予測子・リンク関数

以下の一般化線形モデルで理解するのがわかりやすいです。

g(\mu_i) = \beta_0 + \beta_1x_i

y \sim Normal(\mu_i, \sigma ^2)

\beta_0 + \beta_1x_iが線形予測子で、 g()がリンク関数ですね。本書の例ではリンク関数は恒等関数になってますが、ここを色々弄れるわけですね。

brms

brmsパッケージを使うとstanファイルを書かなくてもお手軽にベイズ統計モデリングが出来ます。ただ、状態空間モデルなど現状未対応なモデルもあるため、stanファイルで記述するベーシックな書き方もちゃんとマスターしといた方が良いようです。

以下に実行イメージを掲載しておきます。*5

library(rstan)
library(brms) # このパッケージ

rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

# データ読み込み
data <- read_csv("hoge.csv")

# モデル作成
result_brms <- brm(
    formula = y ~ X,
    family = gaussian(link = "identity"), # 確率分布
    data = data,
    seed = 1
)

# MCMCサンプルの取得
# as.mcmc(result_brms, combine_chaines = TRUE)

# 事後分布の確認
plot(result_brms)

4部 【応用編】一般化線形混合モデル

ここからは階層ベイズモデルを扱います。3部で扱ったポアソン回帰で 過分散が起きていた場合の対処方法として、一般化線形混合モデルを用いて対処する方法が紹介されています。

固定効果・ランダム効果

3部で推定してきたパラメータ（切片、”天気が釣獲尾数にもたらす効果”、”気温が釣獲尾数にもたらす効果”といった部分）が固定効果です。 一方線形予測子に加えられた調査ごとに変化するようなランダムな影響をランダム効果と言います。

固定効果、ランダム効果両方を考慮するモデルを混合モデルと言います。 つまり一般化線形混合モデルというのは、一般化線形モデルにランダム効果を加味したモデルということですね。

5部 【応用編】状態空間モデル

最後は状態空間モデルです。ローカルレベルからトレンド、周期性、自己回帰モデル、動的一般化線形モデルと様々なモデルを試せるようになっています。

状態モデル・観測モデル

状態空間モデルは”状態モデル”と”観測モデル”から成ります。まず状態モデルは以下のように表せます。

\alpha_t \sim f(\alpha_t|\alpha_{t-1})

つまり、t時点の状態がt-1時点の状態に影響されるということですね。この状態に誤差が加わったものが観測値と考えるわけです。つまり観測モデルは以下のように表せます。

y_t \sim g(y_t|\alpha_{t-1})

ローカルレベルモデル

ローカルレベルモデルは以下のようなものです。

\mu_t = \mu_{t-1} + w_t, w_t \sim Normal(0, \sigma_w ^2)

y_t = \mu_t + v_t, v_t \sim Normal(0, \sigma_v ^2)

w_tを過程誤差、 v_tを観測誤差と言います。 \mu_tは正規ホワイトノイズに従う過程誤差が積み重なっているような系列です。つまりランダムウォークですね。因みにホワイトノイズというのは、期待値が0、分散が一定で、同時刻以外の自己相関が0となっているものです。さらに”正規”なのでこれはi.i.d系列になっているのでした。

基本構造時系列モデル

観測値が「トレンド成分 + 周期成分 + ホワイトノイズ」で表現されるモデルを基本構造時系列モデルと言います。例えば、1次のトレンドに周期kの周期成分が入ったモデルは以下のようになります。

\mu_t \sim Normal(\mu_{t-1}, \sigma_w ^2)

\gamma_t \sim Normal(- \sum_{i=t-(k-1)}^{t-1} \gamma_i, \sigma_s ^2)

\alpha_t = \mu_t + \gamma_t

y_t \sim Normal(\alpha_t, \sigma_v ^2)

\mu_tがトレンド、 \gamma_tが周期成分ですね。周期成分について補足をすると、周期kのk個の周期成分の合計が満たす性質、

\sum_{i=t-(k-1)}^{t} \gamma_i = s_t, s_t \sim Normal(0, \sigma_s ^2)

に注意して、

\gamma_t = - \sum_{i=t-(k-1)}^{t-1} \gamma_i

となります。

全体の感想

今まで統計モデリングの本というと、『みどり本』と『アヒル本』が2強という感じでした。私は『アヒル本』はまだ全部読んではいないんですが、入門という観点だと本書がベストな気がしています。

というのも、とにかく例が豊富で、コード解説も初心者にもわかりやすく工夫されているので、本書を実際に手を動かしながら読み進めることで、stanの基本的な使い方や統計モデリングの勘所が自然と身に付く構成になっていると思います。ただ、当然入門書なので、細かいところ*6は追加で適宜勉強していく必要があります。

とはいえ本書だけでもまあまあお腹一杯感があるんですが（笑）これで満足せず、引き続きベイズ統計モデリングを学んでいければなと思います。

*1:データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学

*2:StanとRでベイズ統計モデリング (Wonderful R

*3:今後この形式が変わる可能性が大いにあります。

*4:こういった面白い関係が出るとどうしても因果関係があるようなストーリーで語りたくなってしまうので。。

*5:このままだと動かないのでデータの読み込みやformulaの部分など適当に変えて下さい。

*6:MCMCサンプリングの細かいアルゴリズムとか

このブログのコンセプト

初めまして、とある分析会社で働いている者です。

学生時代は数学の沼にどっぷり浸かり*1、大学院で金融工学を摘み食いした後、データサイエンティストとして世に放たれました。

当時は、データサイエンティストブーム真っ盛り？で未経験の私にもチャンスを頂けて、それからコツコツとやれることをやってきた結果、何とか首の皮一枚繋がって生き延びることが出来ております。

そんな私ですが、スキルを金融ドメインに振りすぎたり、ウェブアプリ開発に精を出したりとちょいちょい脱線した結果、データサイエンティストとして結構ムラのある人材になってしまったなあと反省している日々です。

そこで、とりあえず今一度原点に戻って、まずはデータサイエンティストとして1.8流くらいを目指そうではないかと思い立ち、せっかくだから ブログにでも記録を残そうと思った次第です。

したがって、当面は修行として取り組んだことや読んだ本のまとめやらをメインコンテンツとして書いていけたらなと思っています。どうぞよろしくお願いします。

*1:尚、出来ない側です