ベイズとMCMCと統計モデルの関係

最終更新:2015年12月2日
最終更新:2016年9月22日

MCMCとは乱数発生アルゴリズムです。ランダムなデータを発生させるアルゴリズムです。
MCMCを使う目的は、統計モデルのパラメタを推定することです。

このページでは、MCMCや、ギブスサンプラー、HMCといったアルゴリズムの詳細には立ち入りません。

この記事では、以下の問いに答えます。
・なぜベイズ統計学に乱数発生アルゴリズムがかかわってくるのか
・なぜ乱数を発生させることで、統計モデルのパラメタが推定できるのか

 

この記事はベイズ推定を応用して状態空間モデルを推定する一連の記事の一つです。
記事の一覧とそのリンクは以下の通りです。
・ベイズ統計学基礎
・ベイズと統計モデルの関係
・ベイズとMCMCと統計モデルの関係
・Stanによるベイズ推定の基礎
・Stanで推定するローカルレベルモデル

スポンサードリンク

 

目次

1．ベイズの定理の復習
2．統計モデルは、ベイズの定理のどこに入ってくるか
3．ベイズでパラメタを推定する
4．確率分布から点推定をする方法
5．EAP推定量を求める
6．積分の代わりとしての乱数発生アルゴリズム

 

1．ベイズの定理の復習

まずは、ベイズの定理と統計モデルについて復習します。
ベイズの定理を知らない方は、こちらをご覧ください。
統計モデルとベイズの関係がわからない方は、こちらをご覧ください。
特にベイズの定理については、わかっていなければ結構つらいので、ぜひ前のページに戻って確認されてください。

ベイズの定理をざっと復習します。
まずは、ベイズ更新の考え方からです。

$$事後確率=事前確率\times{修正項}$$

こちらの式のように、ベイズの定理は、事前確率を事後確率に更新するのでした。
修正項をもう少し丁寧に書くと、以下のようになります。

$$事後確率=事前確率\times{\frac{ある状況で、そのデータが得られる確率}{平均的に、そのデータが得られる確率}}$$

こちらを数式で書くと、以下のようになります。

$$P(θ|X)=P(θ)\times{\frac{P(X|θ)}{P(X)}}$$

ちょっと式変形すると、おなじみのベイズの定理が出てきます。

$$P(θ|X)=\frac{P(X|θ)P(θ)}{P(X)}$$

ただし、θは統計モデルのパラメタであり、Xは得られたデータを表します。
P(θ）は事前確率です。
P(θ|X)は事後確率です。

例えば、以下のような線形回帰モデルを想定します。

ビールの売り上げ ＝ 傾き×気温 + 切片

この時のXには、「気温のデータ」が入ります。
この時のθは「傾き」と「切片」が入ります。
しかし、推定すべきパラメタが2つもあると、少々話が複雑になってしまうので、「切片は３だとわかっている」と仮定しましょう。というわけで、推定すべきパラメタθは傾きとなります。

P(θ）は、例えば「最初は、傾きの値が5~6の間にはいる確率は20％だと考えていた」という事前確率を意味します。
P(θ|X)は、例えば「データXが手に入ったことにより、傾きが5~6の間にはいる確率が40%に変化した」という事後確率を表します。

次に傾きθについて考えます。
θは3になるかもしれないし、4.5かもしれないし、－2.376になるかもしれません。
θは小数点以下の値など様々な値をとりうる変数です。
このような連続型の変数では「θがちょうど1になる確率」も「ちょうど－2.376になる確率」も、どんなものでもすべて確率が0となってしまいます。
というわけで、連続型の変数を扱う場合は、確率の代わりに「確率密度」を使います。

θが入っているときは確率Pを使えないので、ちょっと式を変形します。

$$f(θ|X)=\frac{f(X|θ)f(θ)}{f(X)}$$

Pという記号を使わないようにしただけの修正ですので、あまり気をもまないでください。

2．統計モデルは、ベイズの定理のどこに入ってくるか

次は、線形回帰モデルを、先ほどのベイズの定理のどこに適用するかを見ていきます。
線形回帰モデルは、以下の構造をとるとします。

ビールの売り上げ ＝ θ×気温 + ３

さらに、ビールの売り上げは正規分布に従うと仮定します。
ちょっと強引ですが、その正規分布の分散は4だとわかっていたとします。
正規分布は、期待値と分散の2つのパラメタがわかれば、確率密度がすぐに計算できます。
その確率密度をN(期待値、分散)で表すことにします。

分散は「4」だとわかっています。
期待値は、先ほどの線形回帰モデルから「θ×気温＋３」だと仮定しています。
というわけで、「統計モデルのパラメタがわかっている条件での、データが得られる確率密度」は以下の通りです。

N(θ×気温＋３、４)

これがf(X|θ)の部分になります。
代入すると、こんな感じ。

$$f(θ|X)=\frac{N(θ\times{気温}+3, 4)f(θ)}{f(X)}$$

というわけで、統計モデルを変えると、f(X|θ)の部分が変わってくるということです。
ちなみにf(X|θ)は「尤度」とも呼ばれます。

線形回帰モデルですと、大したことがないのですが、状態空間モデルなんかを対象とする場合は、尤度がとても複雑な式になってしまいます。

 

3．ベイズでパラメタを推定する

ベイズの定理は、事後確率を計算するための公式です。
ベイズの定理を使って計算されるのは「確率（密度）」です。

ベイズ統計学は、統計モデルのパラメタを推定するための道具なのでした。
しかし、直接的に「傾きθは3.5です」と推定できるわけではありません。
ベイズ統計学を用いた場合、推定されるのは「パラメタの事後確率」です。

例えば、先ほどの線形回帰の場合ですと、以下のようになります。
「傾きθが3～3.5の間に入る確率は20％です」
「傾きθが6～6.5の間に入る確率は0.2％です」
これを見れば、傾きは3前後の値をとるのだなと推測がつきますね。
6を超えることは、まずないだろう、とか。

先ほどのベイズの定理を使いますと、推定されるのは事後確率、すなわち、f(θ|X) です。
f(θ|X) とは「データXが手に入ったという条件における、パラメタがθになる確率（密度）」です。

このように、「パラメタが○～×に入る確率が△△％」という推定結果が出ることには、ぜひ注意してください。
推定された結果は「パラメタの確率（密度）」なのです。

また、「パラメタが3~4に入る確率」や「4~5になる確率」「8~9になる確率」など様々な確率を求めることができます。
そうすると、パラメタの確率分布が最終的に得られることになります。
ベイズ推定により得られるものは「パラメタの確率分布」だとご理解ください。

 

4．確率分布から点推定をする方法

パラメタの確率分布が求まるというのは大変結構なことではあるのですが、「傾きの値は○○ですよ」と一言で言いたいときはよくあります。
例えば、予測をしようとしたとき、「θ×気温＋３」でビールの売り上げを計算します。これで「θとは確率分布で……」と言っていては、予測ができません。θの値を一つだけほしいと。

この時には、中央値をとるなど、様々なやり方が使われます。
しかし、一番多いのは、θの期待値をとる方法でしょう。

期待値は「確率×その時の値」の合計値として計算されます。
確率分布さえ分かっていれば、期待値を計算することは難しくありません。

まず、ベイズの定理のおかげで、パラメタの事後確率分布 が計算できています。
パラメタの値はθです。
ということで「確率×その時の値」は 「 f(θ|X) × θ 」となります。

これを合計してやればよいです。

このようにして計算されたθの期待値を「EAP推定量」と呼びます。
（EAPはExpected a Posterioiの略で、事後期待値と訳されます）

 

5．EAP推定量を求める

理論上は、確率分布をベイズの定理を使って計算できたのであれば、そこからパラメタの期待値を計算し、EAP推定量を求めることができるはずです。
できるはずなんですが、これが、実際にやってみるととても難しいのです。

期待値は「確率×その時の値」の合計値です。
しかし、パラメタが小数点以下をとるなど、連続量をとるものであった時には、単純に合計はできません。足し算を無限回しなければならないからですね。
そういう時に使われるのが積分です。
連続量を扱う場合は、シグマをとれないので、積分をするということです。

積分をすると、以下のような計算式になります。

$$EAP推定量　=　\int f(θ|X)\timesθ　dθ$$

積分が苦手な方は「足し算の代わり」という認識で積分記号を眺めてください。
「確率×その時の値」である「 f(θ|X) × θ 」を足し合わせているということです。そしたら期待値が計算できます。

そして、ベイズの定理を使うと、事後確率分布を、ベイズ更新の形で書き表すことができるのでした。

$$EAP推定量　=　\int \frac{f(X|θ)f(θ)}{f(X)}\timesθ　dθ$$

そして、尤度であるところの f(X|θ) は、統計モデルが複雑になればなるほど、ややこしい式になります。

例えば、線形回帰分析の単純な例では以下のようになります。

$$EAP推定量　=　\int \frac{N(θ\times{気温}+3, 4)f(θ)}{f(X)}\timesθ　dθ$$

これくらいであれば、何とかなるのですが、これが状態空間モデルになり、推定すべきパラメタが10個とか20個とかになってくると、もう駄目です。積分の計算ができなくなります。
「積分ができなくなる」というのは、手計算が難しくなるという意味ではありません。
コンピュータを使っても、計算するのが無理になるということです。

 

まとめます。

ベイズの定理を使うと、パラメタの確率分布が求まります。パラメタが直接「2になります」とか「3になります」とかいうのはわかりません。
ベイズの定理とは確率を求める公式であったことをぜひ思い出してください。

とはいえ、パラメタが一つ「点」として求まるよりも、確率分布として求まったほうが便利です。なぜならば、中央値を計算することもできるし、期待値を計算することもできる。様々な統計量を計算することができるからです。
ですが、この計算に問題があります。
式があまりにも複雑になってしまい、積分ができなくなるという問題です。
積分ができなくなると、期待値が計算できません。もちろん、期待値の計算以外にも、様々な弊害があります。普段使っている統計量が軒並み計算できなくなるでしょう。
大変困ります。

そこで、乱数発生アルゴリズムの出番です。

 

6．積分の代わりとしての乱数発生アルゴリズム

積分の代わりに乱数を使います。
ここでのポイントは「確率を数式で求める代わりに、シミュレーションを使って求める」というところです。

例えば「平均4、分散3の正規分布に従う」シミュレーションデータが100万個あったとしましょう。
このデータの平均値はどのようにして計算すればよいでしょうか。
平均4とすでに分かっているじゃないかという突っ込みは置いておきます。
平均を求めたければ、100万個あるシミュレーションデータの平均値をとりますね。

次は「形状母数 が4、尺度母数が2のガンマ分布に従う」シミュレーションデータが100万個あったとします。
形状母数が何かわからない人であっても、ガンマ分布から期待値を計算する公式がわからない人でも、平均値はすぐに計算できますね。
100万個あるシミュレーションデータの平均値をとればよいだけです。

それと同じです。
乱数発生アルゴリズムを使って「パラメタθの事後確率分布に従う」シミュレーションデータを作ります。
例えば、シミュレーションをしてパラメタθを100万個作ったとしましょう。
その100万個あるパラメタθを使えば、「θが3～3.5の間に入る確率」が計算できるし、期待値だって計算できます。
期待値の場合は、何も考えずに、100万個あるθの平均値をとるだけです。こうすればEAP推定量が計算できます。

とはいえ、シミュレーションしてθをたくさん作るというのは簡単ではありません。テキトーにシミュレーションして、事後確率分布と全然違ったものが出てきては困ります。
ですので、MCMCといった少し高度な手法が必要になってきます。

というわけで、積分が難しいベイズの定理の結果を使う場合には、乱数発生アルゴリズムを使ってパラメタを推定するのだよというお話でした。

 

参考文献

基礎からのベイズ統計学: ハミルトニアンモンテカルロ法による実践的入門 何度も紹介していますが、ベイズを学ぶにはやはりこの本がよいと思います。 ベイズ推定の考え方から実践の仕方まで、一通り学ぶことができます。

岩波データサイエンス Vol.1 特集「ベイズ推論とMCMCのフリーソフト」 主にソフトの使い方ですが、著名な方が丁寧に解説してくださっているので、勉強になるところが多いです。 ソフトの使い方以外の部分もぜひ読まれるとよいかと思います。 ベイズ超速習コースという記事がありまして、そちらも参考になります。 読みやすさと内容のバランスが絶妙の本で、お値段も安いので、お勧めです。

リンク

この記事はベイズ推定を応用して状態空間モデルを推定する一連の記事の一つです。
記事の一覧とそのリンクは以下の通りです。
・ベイズ統計学基礎
・ベイズと統計モデルの関係
・ベイズとMCMCと統計モデルの関係
・Stanによるベイズ推定の基礎
・Stanで推定するローカルレベルモデル

スポンサードリンク