ハミルトン系の数値計算 (PRML11章)

PRML (→amazon;パターン認識と機械学習 下) の11章に出てくるハミルトン系の話。ここではハイブリッドモンテカルロのためにハミルトン系を導入して、リープフロッグスキームと呼ばれる数値積分(差分スキーム)を導入している。リープフロッグとかは昔ひと通り勉強したのだけど、久しぶりに見て懐かしかったので思い出しながら/遊びがてら計算してみた。
“ハミルトン系の数値計算 (PRML11章)”の続きを読む

Rからc++コードを呼び出すRcpp+inlineを試す

従来、RをCで拡張するときは「Rの基礎とプログラミング技法(→amazon)」にもあるとおり、R.h をインクルードして SEXP (S EXPression) 型と格闘するイメージだけど(何度か実践で使ったけど気楽にやろう、ってものではなかった)、Rcpp + inline を使うと魔法のように簡単に R と c++ を連携させることができるようだ。

なんとインラインで c++ を書いて、それを直接実行する、というすごい実装(inlineパッケージを使った場合。Rcpp は R.h のラッパー)。

使い道としてはMCMCとかやたら遅いアルゴリズムを実装するときに便利そう(C++の中でRの豊富な乱数生成が使えたりするから)。

他にもループが遅い時にループの深いところを c++ で書き換えてしまうとか。かなり気楽だからこそできるワザ。

勉強を兼ねて適当なテスト関数を書いてみた。
“Rからc++コードを呼び出すRcpp+inlineを試す”の続きを読む

正定値行列と共分散とフィッシャー情報量(と機械学習との関連について少々)


確率統計ででてくる基本的な2つの行列(共分散行列、フィッシャー行列)の半正定値性はほとんど同じ方法で証明できる、ということについて。最後の方は機械学習と関連したすこしマニアックな話題。

関連1 → Cramer-Rao’s theorem
関連2 → 正定値行列と共分散(に関する駄文)

共分散行列の半正定値性

そもそも、共分散行列が半正定値になる理由は、任意のベクトル \(u\) に対して
\begin{align}
u^T\rm{Cov}[\boldsymbol X]u
&=\int u^T(\boldsymbol x-\bar{\boldsymbol x})
(\boldsymbol x-\bar{\boldsymbol x})^TudF(x)\\
&=\int |u\cdot(\boldsymbol x-\bar{\boldsymbol x})|^2dF(x)\ge0
\end{align}
となることから明らか。

フィッシャー情報量の半正定値性

同様に、フィッシャー情報量が半正定値行列になる理由は \(l(x|\theta)=\log f(x|\theta)\) として
\begin{align}
u^TI(\theta)u
&=u^TE[\nabla_\theta l(X|\theta)\cdot\nabla_\theta^Tl(X|\theta)]u\\
&=\int |u\cdot\nabla_\theta l(x|\theta)|^2dF(x|\theta)\ge0
\end{align}
より明らか。よって逆行列が存在すれば(つまりフィッシャー行列が退化していなければ) \(I(\theta)^{-1}\) も正定値。

さらに、任意の行列 \(A\) に対して \(A^TI(\theta)^{-1}A\) も半正定値(\(A\) が退化することもあるので「半」がくっついてくる)。これも上と全く同じ論法で終了。なのでクラメール・ラオのエントリーで出てきた次の式
\begin{align}
\mathrm{Cov}_{\boldsymbol{\theta}}
\left(\boldsymbol{W}(\boldsymbol X)\right)\geq \frac {\partial \boldsymbol{\tau}
\left(\boldsymbol{\theta}\right)} {\partial \boldsymbol{\theta}}
[I\left(\boldsymbol{\theta}\right)]^{-1}
\left( \frac {\partial
\boldsymbol{\tau}\left(\boldsymbol{\theta}\right)}
{\partial \boldsymbol{\theta}}\right)^T
\end{align}
の右辺も左辺も半正定値になることがわかる。

機械学習とフィッシャー情報量

このくらいまで書くとフィッシャー情報行列が退化する、というワクワクする状況についても書きたくなってくる。フィッシャー情報量 \(I(\theta)\) がモデルの真のパラメータのところで退化するケース。ゼロ固有値を持ち、その逆行列は存在しない(形式的には固有値ゼロで割るので逆行列が無限大に発散する)。こういう状況は機械学習でしばしば生じる。

(誤差項に正規分布を仮定して確率モデルとして扱った場合の)多層パーセプトロンとかがこれに該当する。他にはHMMとか隠れノードのあるベイジアンネットとか。hidden state がある場合に退化することが多いっぽい。ちなみに SVM はそもそも確率モデルでないので対象外。

フィッシャー行列が退化すると何がまずいのか、というと、多層パーセプトロンはバックプロパゲーションするとローカルミニマムに落ちるからそもそも最尤推定は難しいのだけれども、直感的に言うと、たとえ最尤推定できたとしても、クラメール・ラオの下限が発散するので、最尤推定量のバラつきも無限大になる、ということを暗に意味している。

これは(悪い意味で)指数型分布族とは著しく異なる性質だ。指数型分布族的なものを正則モデル(regular model)というのに対して、特異モデル(singular model)というらしい。こいつらに対しては AIC などのモデル選択も理論的にはうまくいかない(実用的にはそれなりにうまくいくという話は周りの人たちからはよく聞く)。

そういうわけで、点推定はうまくいかないので、ベイズ推定をすることになる。ベイズ推定なら特異モデルに対しても漸近的な性質を導くことができる。くわしいことは特異モデルの機械学習のベイズ推定の理論の本 “Algebraic Geometry and Statistical Learning Theory” に書いてある(難易度:個人的な感想としては激高)。

正定値行列と共分散(に関する駄文)

昨日、クラメール・ラオの下限について書いたが、そもそもなんで
\begin{align}
\mathrm{Cov}_{\boldsymbol{\theta}}
\left(\boldsymbol{W}(\boldsymbol X)\right)\geq \frac {\partial \boldsymbol{\tau}
\left(\boldsymbol{\theta}\right)} {\partial \boldsymbol{\theta}}
[I\left(\boldsymbol{\theta}\right)]^{-1}
\left( \frac {\partial
\boldsymbol{\tau}\left(\boldsymbol{\theta}\right)}
{\partial \boldsymbol{\theta}}\right)^T
\end{align}
だと共分散行列の下限を定めたことになるの?という話について。

一次元の場合は
\begin{align}
V[X]\ge V[Y]
\end{align}
ならば X のバラつきは Y より大きい、ということはわかるけど、
\begin{align}
\rm{Cov}[\boldsymbol X]\ge A
\end{align}
という場合、つまり \(\rm{Cov}[\boldsymbol X]-A\) が(半)正定値行列のとき、この式から \(\rm{Cov}[\boldsymbol X]\) について何が言えるのだろうか?

以下、答え。\(\rm{Cov}[\boldsymbol X]\) の対角成分に注目する。これはそれぞれ確率変数ベクトル \(\boldsymbol X\) の第 i 成分の分散に相当する。\(\rm{Cov}[\boldsymbol X]-A\) が正定値行列という仮定より、この行列の対角成分は0以上になる(後述)。つまり \(V[X_i]\ge a_{ii}\) がいえる。なので最初に戻って
\begin{align}
\mathrm{Cov}_{\boldsymbol{\theta}}
\left(\boldsymbol{W}(\boldsymbol X)\right)\geq \frac {\partial \boldsymbol{\tau}
\left(\boldsymbol{\theta}\right)} {\partial \boldsymbol{\theta}}
[I\left(\boldsymbol{\theta}\right)]^{-1}
\left( \frac {\partial
\boldsymbol{\tau}\left(\boldsymbol{\theta}\right)}
{\partial \boldsymbol{\theta}}\right)^T
\end{align}
から言えることは、

  • \(\boldsymbol W(\boldsymbol X)\) の成分ごとの分散は右辺の行列の対角成分の値以上になる

となってやっぱり下限を定めていることになる。というあたりまえのことをなんとなくやる気なさげに書いてみました。おしまい。

おまけ

半正定値行列の任意の対角成分は0以上になることは、半正定値性の定義、つまり任意のベクトル u に対して
\begin{align}
 u^TAu\ge0
\end{align}
となることを考えればすぐわかる。つまり u を第i成分だけ1であとは0みたいなベクトルにすれば対角成分だけ取り出されて \(a_{ii}\ge0\) がわかる。

Cramer-Rao’s theorem

Cramer-Rao’s theorem について.簡単に言うと確率モデルのパラメータを推定する際の精度の理論的限界を定める定理(もちろんこの定理を適用するためにはいくつかの仮定=前提条件があるのですが).1パラメータのケースは 『Casella & Berger 本』 に書いてあるけれど,パラメータが1次元より大きなときについては書いてなかったので考えてみた.

(自分のメモをブログに転載しただけなのであまりフレンドリーな書き方ではないかもしれませんが,いつか誰かの役に立つかもしれないので残しておきます.昔に書いた関連エントリはコチラ

考えている状況とか、記号の定義とか

\(\boldsymbol X\) を \(f(\boldsymbol X|\boldsymbol\theta)\) に従う確率変数とする.\(\boldsymbol{\theta}\) は分布の母数.\(\boldsymbol\tau(\boldsymbol\theta)\) は母数を一対一変換するような微分可能な関数.\(\boldsymbol W(\boldsymbol X)\) は \(\boldsymbol\tau(\boldsymbol\theta)\) の不偏推定量,つまり
\begin{align}
E_{X|\theta}(\boldsymbol W(\boldsymbol X))&=\boldsymbol\tau(\boldsymbol\theta)
\end{align}
が成立するとする.ようするに \(\boldsymbol W\) は観測されたデータをモデルパラメータに変換する写像(アルゴリズム)と考えればよい.\(\boldsymbol W(\boldsymbol X)\) は確率変数なので,バラつきがある.つまり推定されたモデルパラメータ自身が確率変数になっている.このバラつきが大きいと精度の良い推定は望むべくも無い.そこでパラメータのバラつきを評価するの理論的な限界(下限)を定めるための定理が Cramer-Rao.推定量といってもいろいろあって,「データがどんな値をとってもパラメータは0!」と決め打ちしてしまうアルゴリズムも考えられる.こういう推定量はバラつきはゼロだけど精度がいいとは言えない.不偏推定量という制約条件を課しているのはこのようなおかしな推定量を排除して「筋の良いもの」だけを相手にするための仕掛けと考えればいい.

この定理があってはじめて,「よい推定方法とは何か?」という問いに答えられるようになる.すなわちその答えは「不偏でかつバラつきの理論的限界を達成するような推定量」こそが最良である,となる.

定理;一般化されたクラメール・ラオの下限

\(I(\boldsymbol\theta)\) をフィッシャー情報量とするとき,
\begin{align}
\mathrm{Cov}_{\boldsymbol{\theta}}
\left(\boldsymbol{W}(\boldsymbol X)\right)\geq \frac {\partial \boldsymbol{\tau}
\left(\boldsymbol{\theta}\right)} {\partial \boldsymbol{\theta}}
[I\left(\boldsymbol{\theta}\right)]^{-1}
\left( \frac {\partial
\boldsymbol{\tau}\left(\boldsymbol{\theta}\right)}
{\partial \boldsymbol{\theta}}\right)^T
\end{align}
が成り立つ.ただし,\(A\ge B\) は \(A-B\) が半正定値行列であることを意味している.

以下、証明っぽいもの。 “Cramer-Rao’s theorem”の続きを読む

中国のアキバ(中関村)の風景

閑話休題。中国のアキバといわれる中関村の風景を撮りました。

中関村駅は4号線の北京大学の隣の駅。北京の北西部には大学が集中している。

かなりでかいビルに電子機器の小売店がつまっている。ラジオ会館みたいな感じ。

“中国のアキバ(中関村)の風景”の続きを読む