ベイズ統計の歴史とSTANを使ったベイズ推定

ベイズ統計の歴史とSTANを使ったベイズ推定

前回までに作成した、サーバー/DBのプラットフォームに機械学習の結果を組み合わせると人工知能系のアプリケーションができる。今回はそれら機械学習の中でも、モデリング系に分類されるベイズ推定について述べたいと思う。

まず最初の参考図書が、「異端の統計学ベイズ」。Sharon Bertsch McGrayneによる「The Theory That Would Not Die: How Bayes’ Rule Cracked the Enigma Code, Hunted Down Russian Submarines, & Emerged Triumphant from Two Centuries of Centuries of Controversy」の邦訳で、ベイズの法則が以下にして生まれ、どのように発展してきたかの歴史を紐解いたものとなる。

ベイズ統計の歴史は、1740年代のトーマス・ベイズから始まる。神学と数学を学んでいたベイズは、当時の「あらゆるものの第一原因は神である」という思想に基づいて、因果の根本にある神=宇宙に潜む秩序の存在を数学的に明らかにしようとした。

ベイズより少し前の時代に複素数と三角関数の定理であるド・モアブルの法則を見出した数学者アブラーム・ド・モアブルが原因から結果へと考えを進めることで確率に関する問題を解いたのに対して、ベイズは逆の方向である結果から原因への逆確率の問題を解こうとした。正に見えているものから、その根本原因(宇宙の秩序=神)を導き出そうとしたのである。

それらを解くための手段として、出発点として、とりあえず何らかの数値(ベイズのいうところの「推測値」)をでっち上げ(主観的確率とも呼ばれる)、情報が得られた時点でその数値を修正するという解決法を思いつき、具体的な思考実験として、後ろにあって見ることの出来ない四角いテーブルにボールを投げ、最初に落ちた場所に対して、その後投げたボールが最初のボールよりも右に行ったか左に行ったかの情報を得ていく事で推定の精度をあげていくことを考えた。

このようなケースの場合、ボールを投げる数を増やすことにより、新たに得られる情報が増えて、最初に投げたボールが落ちたと思われる場所の推定範囲が狭まっていく。これは概念的には、客観的な情報に基づいて自分の意見を変えると言う極々当たり前の話でもある。

ここで情報と確率を結びつける。シャノンによる情報理論では、情報の量は不確かさが多い(確率が小さい)ものが大きくなり、少ないものが小さくなると定義している。別の言い方をすると、これから起こることが100%分かっている情報を受け取ってもありがた味は感じないが(価値はない)、起こる確率がの小さい事象についての情報を受け取る(それが起きることが分かっていれば)ありがたさがを感じる(価値がある)と言うことにになる。つまり情報量と確率は関連すると言うことになる。

この考え方を進めて、二つの情報を足し合わせる場合を考える。 例えば1/2と1/4でそれぞれ起こる確率の情報があるとした場合、これら2つの情報を足し合わせてると言うことは、具体的な例だと4枚のカードから抜いたものが偶数である情報は確率1/2で、その数字が特定のもの例えば2とか3とかである情報は1/4、それらを足し合わせたものは特定の数字でかつ偶数なので2とか4である情報つまり1/2×1/4=1/8の確率である情報だと言える。つまりより細かで詳細な情報は様々な情報の足し算であり、確率的にはそれらを掛け算した希少な確率の情報となる。これは情報の足し算が確率の掛け算となることを示している。

以上よりベイズの情報を集めていくという概念は、それぞれの起こる確率の掛け算で表されることになる(前提として集める情報の確率は計算できるものとなる)。

情報と確率を結びつけることができたので、更に結果と原因の概念を確率で表す。それらの概念はある原因の元で結果が起こると言い換えることができる。ここで条件付き確率を、事象Aが起こったと言う条件の元で事象Bが起こると言う確率と定義して、 式として

\[P(B|A)=\frac{P(A\cap B)}{P(A)}\]

と定義すると、Aを原因、Bを結果としてこの式で表すことができる。ここで逆の方向から見たP(A|B)=P(B∩A)/P(B)を考え、確率の積の法則よりP(A∩B)=P(B∩A)であることから

\[P(B|A)P(A)=P(A\cap B)=P(B\cap A)=P(A|B)P(B)\]

と表される。ここで原因の確率Aは分からないので、他の確率で表すと以下のようになる。

\[P(A)=P(A|B)P(B)+P(A|\overline B)P(\overline B)\]

これらを用いてベイズの式として以下が導出される。

\[P(B|A)=\frac{P(A|B)P(B)}{P(A|B)P(B)+P(A|\overline B)P(\overline B)}\]

この「異端の統計学ベイズ」に記載されているベイズの逐次的なアプローチとしては、チューリングがエニグマの解読の為や、海に消えた潜水艦や水爆を探したり、あるいは意思決定や、保険数理、がんの要因の発見等がある。これらはあくまでも手作業あるいはそれらを効率化したものでしかなく、少ないパラメータの世界でしか利用できなかった。

これに対して、コンピューターを用いたアプローチが1980年代より行われ始める。多数のパラメータのあるモデルへの対応として階層ベイズやギブスサンプラー等のマルコフ連載モンテカルロ法 (MCMC) 等を使う事で一度に相手にするパラメータを限定して逐次計算していくと言う手法が取られ始める。また機械学習のキモの一つである関数の積分による最適化の手法を組み合わせることで、大きな飛躍が得られた。実際に使われた領域としては、故障分析や病気の要因の抽出、そしてデジタル画像修復等の画像処理の分野等がある。

このMCMCをコード化したものが、BUGS(Bayesian inference Using Gibbs Sampling)やJAGS(Just Another Gibbs Sampler)、Stan(RStan)、PyMC等のオープンソースソフトウェアとなる。参考図書としては、BUGSでは「ベイズ統計で実践モデリング: 認知モデルのトレーニング」

JAGS等では「ベイズ推論とMCMCのフリーソフト」、

「ベイズ統計モデリング―R,JAGS,Stanによるチュートリアル― 原著第2版」、

Stanに関しては「StanとRでベイズ統計モデリング」、

「実践Ｄａｔａ　Ｓｃｉｅｎｃｅシリーズ　ＲとＳｔａｎではじめる　ベイズ統計モデリングによるデータ分析入門」、

またベイズは上記のように統計のアプローチの参考書が多いが、機械学習に特化したものとして「ベイズ推論による機械学習入門」等がある。

機械学習のアプローチとしては、ポアソンやガウス分布にモデルを仮定してそのパラメータを求める混合モデル(mixture model)近似推論、線形次元削(linear dimensionality reduction)や、”非負値行列因子分解（NMF）の概要とアルゴリズム及び実装例について“で述べている非負値行列因子分解(Nonnegative matrix factorization NMF)、隠れマルコフモデル(Hidden Markov model HMM)、トピックモデル、テンソル分解(tensor factorization)、ロジスティック回帰(logistic regression)、ニューラルネットワーク等がある。

実際のツールを使ったものとしてはStanを使ったアプローチがある。StanはMCMCを使ってサンプルデータを生成する、いわば乱数的なシミュレーションの一種である。通常の機械学習は主にモデルと実データのギャップを最小化するために関数を最適化していくのに対して、モデリング系は最初に設定したモデルのパラメータをシミュレーションしながら最適化していく為、学習データが少ないケースで強みを発揮する。

実際にStanを使うには、Rの’RSatn’と言うツールをインストールして利用する。

install.packages('rstan', dependencies=true)

上記はmodelとして、Y[n]が平均a + bX[n]・標準偏差sigmaの正規分布から生成されたものであることを宣言している。これを実行する為には以下のようなコードを書く。

library(rstan)

d <- read.csv(file='input/data.txt')
data <- list(N=nrow(d), X-d$X, Y=d$Y)
fit <- stan(file='model/stan-file.stan', data=data, seed=1234)

stanのコードは以下のようになる。

data {
データYの宣言
}
parameters {
サンプリングしたいパラメータθの宣言
}
modle {
尤度p(Y|θ)の記述
事前分布p(θ)の分布
}

model部は様々な確率分布で構成される。例えば、ある商品を買うか買わないか、薬が効くかどうかなどの結果が二値で表される場合に適用されるベルヌーイ分布や二項分布、品目がN種類あり、それぞれの品目の故障率の分布を考える場合に適用されるベータ分布、文書の中のトピックの分布を求めるのに使われるカテゴリカル分布やディリクレ分布等々。

data部はmodel部で使う、データの個数とデータの本体Yの宣言。parameters部は推定すべきモデルの中のパラメータを宣言する。確率変数であるので分散や平均、標準偏差になる。

最もシンプルな単回帰モデルでは以下のようになる。

data {
  int N ;
  real X[N];
  real Y[N];
}
parameters {
   real a;
   real b;
   real<lower=0> sigma;
}
model {
  for (n in 1:N) {
     Y[n] = normal(s + b*x[n], sigma):
}}

上記はmodelとして、Y[n]が平均a + bX[n]・標準偏差sigmaの正規分布から生成されたものであることを宣言している。これを実行する為には以下のようなコードを書く。

library(rstan)

d <- read.csv(file='input/data.txt')
data <- list(N=nrow(d), X-d$X, Y=d$Y)
fit <- stan(file='model/stan-file.stan', data=data, seed=1234)