乱数発生アルゴリズムによるベイズ推定

モデルの推定

ベイズの役割

ベイズはあくまで、統計モデルを推定する道具のひとつ。例えば「ベイズという名前の統計モデルがある」という認識は誤り。

「統計モデルを推定する」という行為

「統計モデルを推定する」という行為は、以下の2つの作業に分けられる。

統計モデルの型を決める
統計モデルのパラメータを推定する

統計モデルの推定(1) 統計モデルの型を決める

統計モデルの「型」とは、正確な用語ではないのが、その真意はわかるんじゃないのかなと。例えば、線形回帰モデルという統計モデルの「型」がある。さらに、このモデルは「気温の上昇に伴い、ビールの売り上げが線形に変化する」というような構造を持っているとする。この構造のことをいう。このモデルを数式っぽく書くと、以下のようになる。

ビールの売り上げ $=$ 傾き $\times$ 気温 $+$ 切片

この式の場合、気温が1度上昇すると、ビールの売り上げは傾きの値分だけ増加する。このように、「興味のある値(ここでいうビールの売り上げ)はどのような構造で変化するか」を決定することが、統計モデルの「型」を決定することになる。

統計モデルの推定(2) 統計モデルのパラメータを推定する

気温が1度上昇すると、ビールの売り上げは傾きの値分だけ増加することがわかったとする。すなわち、モデルの「型」(モデルを表す式)が決定しているということになる。次に気になるのはその「傾き」と「切片」の値。こういった「傾き」や「切片」といった値が、いくらになるのかを推定する行為を「パラメータの推定」と呼ぶ。

ベイズが関わるのは、この部分で、ベイズは「統計モデルのパラメータを推定するための道具」であるということ。逆に言えば、ベイズ統計学は、統計モデルの「型」を決める際にはあまり登場しない(もちろん状況によるが)。

なぜ統計モデルの推定は2段階に分かれているのか

統計モデルを推定するソフトなりライブラリなりパッケージなりがあったとする。そこにデータを入れれば、全自動で統計モデルが出来上がる…ということは、少なくとも現在、ありえない。それは、統計モデルのパラメータを推定するのは、コンピュータの仕事であるが、統計モデルの「型」を決めるのは、人間の仕事であるから。

したがって、コンピュータは統計モデルの「型」を自動で決めることができないため、コンピュータだけで、全自動で統計モデルを作成することは、現状ではできない。

例えば、ビールの売り上げ変動の構造に、どれくらいのパターンがあるだろうか。データが与えられたとき、我々の目の前に現れるのは、無限とも可能性があるはずだ。その無限のバリエーションの中から一つを選ぶ作業こそが、モデルを推定するという行為。 (私個人としては、この作業が一番楽しいところだと思う。)

コンピュータができるのは、「選択肢が与えられた状態で最善を選ぶこと」のみ。そして、その選択肢を与えるのは、人間。コンピュータに任せられるところ、人間が決めなければならないところ、これらを別々に作業しなくてはならないので、統計モデルを推定する際には、2段階のステップを踏む必要がある。

なぜ状態空間モデルの推定に、ベイズが使われるのか

モデル推定において、人間は、現象に対して、人間が理解できる「型」あるいは「構造」を決めなければならない。そして、状態空間モデルは、統計モデルの「型」の一種。

状態空間モデルは、非常に多くの種類のデータに適用することができる、柔軟なモデルの「型」で、ある。人間が理解しやすい構造で、直感に合うモデルを作ることのできる「型」である。特に、複雑な現象を扱う場合に、状態空間モデルはその威力を発揮する。

しかし、このような幅広く使えるモデルの場合、パラメータを推定するのが難しい。そこで、ベイズの出番です。ベイズ統計学、そしてMCMCを使うと、今まで困難だった複雑なモデルであっても、パラメータを(ある程度)正しく推定することができる。ベイズ統計は、便利な道具なのだ。オイシイネ。

ベイズと乱数発生アルゴリズム

ベイズの定理

\begin{aligned} P\left(\theta|X\right) &= \frac{P\left(X|\theta\right)P(\theta)}{P(X)} \end{aligned}

例題

例えば、ビールの売り上げを以下のような線形回帰モデルで表現したとする。

ビールの売り上げ = 傾き\times 気温 + 切片

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

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

傾き $\theta$ は、連続変数をとるため、どんな $\theta$ の値でもすべて確率が0となってしまう。そこで、確率の代わりに「確率密度」を用いる。 $\theta$ が入っているときには、確率 $P$ を使えないので、以下のようにちょっと式を変形する (ただ、Pという記号を使わないようにしただけの修正)。

\begin{aligned} f\left(\theta|X\right) &= \frac{f\left(X|\theta\right)f(\theta)}{f(X)} \end{aligned}

統計モデルへのベイズの定理の適用

上記の例題によれば、線形回帰モデルは、以下の構造をとる。

ビールの売り上げ = \theta\times 気温 + 3

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

分散は4だとわかっていて、ビールの売り上げの期待値は、上記のように線形回帰モデルを仮定している。したがって、「統計モデルのパラメータがわかっている条件の下での、データが得られる確率密度」は以下のようになる。この $f(X|\theta)$ を尤度と呼ぶ。

\begin{aligned} f\left(気温|\theta\right) &= \mathcal{N}(\theta\times 気温 + 3, 4) \end{aligned}

$気温=X$ なので、ベイズの定理に代入すると、以下のようになる。

\begin{aligned} f\left(\theta|気温\right) &= \frac{\mathcal{N}(\theta\times 気温 + 3, 4)f(\theta)}{f(気温)} \end{aligned}

ベイズでのパラメータ推定

ベイズの定理は、事後確率を計算するための式なので、得られる値は、確率 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)$ は連続関数なので、総和は積分で表す。

\begin{aligned} E\left(\theta\right) &= \int f\left(\theta|気温\right)\theta d\theta \end{aligned}

このようにして計算された $\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}

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

\begin{aligned} E\left(\theta\right) &= \int \frac{\mathcal{N}(\theta\times 気温 + 3, 4)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手法は、ベイズ統計という統計学の「考え方」を実際のデータ解析に応用する際に使われる、アルゴリズムの一種。