乱数発生アルゴリズムによるベイズ推定¶
モデルの推定¶
ベイズの役割¶
ベイズはあくまで、統計モデルを推定する道具のひとつ。例えば「ベイズという名前の統計モデルがある」という認識は誤り。
「統計モデルを推定する」という行為¶
「統計モデルを推定する」という行為は、以下の2つの作業に分けられる。 1. 統計モデルの型を決める 1. 統計モデルのパラメータを推定する
統計モデルの推定(1) 統計モデルの型を決める¶
統計モデルの「型」とは、正確な用語ではないのが、その真意はわかるんじゃないのかなと。 例えば、線形回帰モデルという統計モデルの「型」がある。 さらに、このモデルは「気温の上昇に伴い、ビールの売り上げが線形に変化する」というような構造を持っているとする。この構造のことをいう。 このモデルを数式っぽく書くと、以下のようになる。
ビールの売り上げ \(=\) 傾き\(\times\) 気温 \(+\) 切片
この式の場合、気温が1度上昇すると、ビールの売り上げは傾きの値分だけ増加する。 このように、「興味のある値(ここでいうビールの売り上げ)はどのような構造で変化するか」を決定することが、統計モデルの「型」を決定することになる。
統計モデルの推定(2) 統計モデルのパラメータを推定する¶
気温が1度上昇すると、ビールの売り上げは傾きの値分だけ増加することがわかったとする。 すなわち、モデルの「型」(モデルを表す式)が決定しているということになる。 次に気になるのはその「傾き」と「切片」の値。 こういった「傾き」や「切片」といった値が、いくらになるのかを推定する行為を「パラメータの推定」と呼ぶ。
ベイズが関わるのは、この部分で、ベイズは「統計モデルのパラメータを推定するための道具」であるということ。 逆に言えば、ベイズ統計学は、統計モデルの「型」を決める際にはあまり登場しない(もちろん状況によるが)。
なぜ統計モデルの推定は2段階に分かれているのか¶
統計モデルを推定するソフトなりライブラリなりパッケージなりがあったとする。 そこにデータを入れれば、全自動で統計モデルが出来上がる…ということは、少なくとも現在、ありえない。 それは、統計モデルのパラメータを推定するのは、コンピュータの仕事であるが、 統計モデルの「型」を決めるのは、人間の仕事であるから。
したがって、コンピュータは統計モデルの「型」を自動で決めることができないため、 コンピュータだけで、全自動で統計モデルを作成することは、現状ではできない。
例えば、ビールの売り上げ変動の構造に、どれくらいのパターンがあるだろうか。 データが与えられたとき、我々の目の前に現れるのは、無限とも可能性があるはずだ。 その無限のバリエーションの中から一つを選ぶ作業こそが、モデルを推定するという行為。 (私個人としては、この作業が一番楽しいところだと思う。)
コンピュータができるのは、「選択肢が与えられた状態で最善を選ぶこと」のみ。 そして、その選択肢を与えるのは、人間。 コンピュータに任せられるところ、人間が決めなければならないところ、これらを別々に作業しなくてはならないので、 統計モデルを推定する際には、2段階のステップを踏む必要がある。
なぜ状態空間モデルの推定に、ベイズが使われるのか¶
モデル推定において、人間は、現象に対して、人間が理解できる「型」あるいは「構造」を決めなければならない。 そして、状態空間モデルは、統計モデルの「型」の一種。
状態空間モデルは、非常に多くの種類のデータに適用することができる、柔軟なモデルの「型」で、ある。 人間が理解しやすい構造で、直感に合うモデルを作ることのできる「型」である。 特に、複雑な現象を扱う場合に、状態空間モデルはその威力を発揮する。
しかし、このような幅広く使えるモデルの場合、パラメータを推定するのが難しい。 そこで、ベイズの出番です。 ベイズ統計学、そしてMCMCを使うと、今まで困難だった複雑なモデルであっても、パラメータを(ある程度)正しく推定することができる。 ベイズ統計は、便利な道具なのだ。オイシイネ。
ベイズと乱数発生アルゴリズム¶
ベイズの定理¶
例題¶
例えば、ビールの売り上げを以下のような線形回帰モデルで表現したとする。
ビールの売り上げ = 傾き\times 気温 + 切片
この時のベイズの定理の\(X\)には、「気温のデータ」が入り、 \(\theta\) は「傾き」と「切片」が入る。 しかし、推定すべきパラメータが2つもあると、少々話が複雑になってしまうので、 「切片は3だとわかっている」と仮定する。というわけで、推定すべきパラメータ \(\theta\) は傾きとなる。
\(P(\theta)\)は、例えば最初は、「傾きの値が5から6の間にはいる確率は20％だ」という事前確率を意味する。 \(P(\theta|X)\)は、例えばデータXが手に入ったことにより、「傾きが5から6の間にはいる確率が40%」に変化したという事後確率を表す。
傾き\(\theta\)は、連続変数をとるため、どんな \(\theta\)の値でもすべて確率が0となってしまう。 そこで、確率の代わりに「確率密度」を用いる。 \(\theta\)が入っているときには、確率\(P\)を使えないので、以下のようにちょっと式を変形する (ただ、Pという記号を使わないようにしただけの修正)。
統計モデルへのベイズの定理の適用¶
上記の例題によれば、線形回帰モデルは、以下の構造をとる。
ビールの売り上げ = \theta\times 気温 + 3
さらに、ビールの売り上げは、正規分布に従うと仮定する。 少々強引だが、さらにその正規分布の分散は4だとわかっていると仮定する。 正規分布は、期待値と分散の2つのパラメータがわかれば、確率密度がすぐに計算できるので、 その確率密度を\(\mathcal{N}(期待値, 分散)\)で表すこととする。
分散は4だとわかっていて、ビールの売り上げの期待値は、上記のように線形回帰モデルを仮定している。 したがって、「統計モデルのパラメータがわかっている条件の下での、データが得られる確率密度」は以下のようになる。 この \(f(X|\theta)\)を尤度と呼ぶ。
\(気温=X\)なので、ベイズの定理に代入すると、以下のようになる。
ベイズでのパラメータ推定¶
ベイズの定理は、事後確率を計算するための式なので、得られる値は、 確率 or 確率密度 。 ベイズ統計学は、統計モデルのパラメータを推定するための道具ではあるが、直接的に「傾き\(\theta\)は3.5です」などと推定できるわけではない。
上記の線形回帰の場合だと、以下のようになる。
| 傾き \(\theta\) がとる区間 | 傾き \(\theta\) が区間に入る確率 |
|---|---|
| \(2.5 - 3.0\) | \(15\%\) |
| \(3.0 - 3.5\) | \(20\%\) |
| \(3.5 - 4.0\) | \(10\%\) |
| \(4.0 - 4.5\) | \(5\%\) |
| \(4.5 - 5.0\) | \(1\%\) |
これを見ると、傾きは3.0前後の値だろう、とか、5.0を超えることはまずないだろう、などと推測がつく。 このように、ベイズの定理を使うと、「データ\(X\)が手に入ったという条件における、パラメータが\(\theta\)になる確率密度」を推定できる、 すなわち、事後確率 \(f\left(\theta|X\right)\)を推定できる。
また、様々なパラメータの値に対する確率密度を推定することができるので、 最終的には、ベイズ推定によりパラメータの確率密度分布を得ることができる。
EAP推定量¶
確率密度分布からの最適パラメータ推定する。 ベイズ推定により、パラメータの確率密度分布が求まるが、ときに、最適なパラメータの値はいくらと一言で言いたいことはありがち。 例えば、予測をしようとしたとき、ビールの売り上げは、 \(\theta\times気温+3\)により計算できるが、 \(\theta\)は確率密度分布なので… となって予測ができない。 \(\theta\)の値を一意に決定したい。ではどうするか。
この時、中央値をとるなど、様々なやり方が使われるが、一般に多いのは、 \(\theta\)の期待値をとる方法である。 期待値は「確率×その時の値」の総和として計算される。 確率密度分布さえ分かっていれば、期待値を計算することは容易。
まず、ベイズの定理により、パラメータの事後確率分布が計算できているので、 パラメータ\(\theta\) 期待値は以下のようにして求まる。ただし、 \(f(\theta)\) は連続関数なので、総和は積分で表す。
このようにして計算された\(\theta\)の期待値 \(E(\theta)\)を「EAP推定量」と呼ぶ。 (EAP : Expected a Posterioriの略で、事後期待値と訳される。)
上記の式をベイズ更新の形で書き換える。 $$ \begin{aligned} E\left(\theta\right) &= \int \frac{f\left(X|\theta\right)f(\theta)\theta}{f(X)} d\theta \end{aligned} $$
例えば、線形回帰分析の単純な例では以下のようになる。
このように単純な式で表せる尤度であれば、この積分が困難になることはないが、 より複雑な、例えば状態空間モデルなどになると、期待値を軒並み計算できなくなる。
乱数発生アルゴリズム¶
積分の代わりに乱数発生アルゴリズムを用いる。 ここでのポイントは、 「確率を数式で求める代わりに、シミュレーションを用いて求める」 というところ。
例えば「平均4、分散3の正規分布に従う」シミュレーションデータが100万個あったとする。 このデータの平均値はどのようにして計算すればよいか(平均4とすでに分かっているじゃないかとツッコミたくても口をつぐめ)。 平均を求めたければ、100万個あるシミュレーションデータの平均値をとればよい。
次は「形状母数 が4、尺度母数が2のガンマ分布に従う」シミュレーションデータが100万個あったとする。 形状母数が何かわからない人であっても、ガンマ分布から期待値を計算する公式がわからない人でも、 100万個あるシミュレーションデータの平均値をとれば、平均値はすぐに計算できる。
これらと同じように考えればよい。 乱数発生アルゴリズムを使って「得られた事後確率分布に従うパラメータ \(\theta\) 」のシミュレーションデータを100万個作成する。 その100万個あるパラメータ \(\theta\) を使えば、 \(\theta\)が\(3-3.5\)の間に入る確率や期待値計算できる。 期待値の場合は、何も考えずに、100万個ある\(\theta\)の平均値をとるだけ。こうすればEAP推定量が計算できる。
とはいえ、シミュレーションして \(\theta\) をたくさん作るというのは簡単ではない。 むやみにシミュレーションして、事後確率分布と全然違ったものが出てきては困る。 なので、Markov Chain Monte-Carlo (MCMC) 手法といった少し高度な手法が必要になる。 MCMC手法は、ベイズ統計という統計学の「考え方」を実際のデータ解析に応用する際に使われる、アルゴリズムの一種。