こんにちは、SI部の杉です。

約一年ぶりの統計解析に関する投稿となります。
Rを用いて色々な解析手法をご覧に入れました。
今回は、ベイズ統計による解析方法を紹介します。

ベイズ統計だけでかなり広い分野であるため、なるべく端的にわかりやすく説明したいと思います。
まずは例によってベイズ統計とは何かを踏まえたうえで実際の解析方法に移っていきます。

ベイズ統計とは何か?の前に統計学のおさらい

統計学には大雑把に分けて、

  1. 頻度主義
  2. ベイズ主義

という概念に分けられます。

1. 頻度主義 ・・・
大まかに表現するとランダムな事象が生起・発生する頻度をもって『確率』と定義する考え方です。
言い換えると、母集団中の割合や確率分布で決まる観測されたデータが『繰り返して起き続ける』という考えです。

2. ベイズ主義 ・・・
確率を主観的なもの、個人の信念の度合いとして定義するという考えです。
言い換えると、結果が確率を決めるということになります。

コインを投げて表裏が出る確率は、1/2 ですが、これは従来の頻度主義によるものです。
真っ平らなコインであれば、表が何回か続けて出ても、表裏の出る頻度は同じになります。
(だって表と裏しかないから!)

対して、例えばレストランのボーイがお昼の来客をカウントしていたとしましょう。
そのレストランには男性客が4人続けて来店してきたとします。
次に来る客は男性か女性か?と考えるときに、このボーイがベイジアン(ベイズ主義者)であれば、
今まで男性が来ているから次も男性が来る確率の方が高いと考えるでしょう。
ベイズ統計をひと言で表現してしまうと、『その時点での情報を元にしたその時点での確率』であるとも言えます。

もともとベイズ統計とは、統計学者であるトーマス・ベイズが示した『ベイズの定理』を統計建てた理論のことです。
このベイズの定理を利用して推定(ベイズ推定という)を行います。
※推定に関しては次章で紹介します。

hindo_shugi bayes_shugi

ベイズ統計解析は、主観を用いて確率を導くということになります。
そして、この主観とは確率分布のことです。

ベイズの定理

ベイズの定理について説明します。
先ほどのレストランの来客が男性である確率を、P(男|来客者)と表記します。
(来客者がいたとして、その人が男である確率です)

ベイズの定理とは、上記の条件付き確率を求めるための計算式です。

bayes_theorem

分母のP(A)とは、Aが発生する確率(上記の例だと来客者が来る確率となります)。
分子のP(A|B)とは、Bの条件のもとAが発生する確率(上記例だと男性の来客者が来る確率となります)。
P(B)は男性である確率ですね。
また、ベイズ推定においてベイズの定理を利用すると、以下の関係が導き出されます。

  • パラメーターθを用いて定義される確率モデルf(y|θ)を考える。
    ※yは実際の観測データ
  • 事前分布π(θ)を与えた時に以下の関係が成り立つ。

bayes_theorem_2

要するに、ベイズ推定とはモデリングを行うことにより上記の等式(特に左辺の事後分布)を導くということになります。
今回は、マルコフ連鎖モンテカルロ法とStanを用いてモデリングを行います。

Rを用いてベイズ統計解析を行うにあたって

モデリングとパラメータの推定に関して、計算の比重が高いため、Rで解析を行うことは理にかなっていると言えます。

  • モデリング ・・・ 観測データをうまく説明できるようにパラメーターの値を決めてあげること
    ※言い換えれば観測データからノイズを取り除き、一定の法則を見つけ出して抽象化すること

   ⇒ 以前取り上げた回帰分析も、プロットから右肩上がりの直線に近似しました。
これも大雑把に言うとモデリングの一種と言えます。

  • パラメータの推定 ・・・ 当てはまりの良さを評価すること(このデータはどれだけ信頼できるか?など)
    ⇒ 以前取り上げた回帰分析の直線近似がどれだけ当てはまっているかを指標するものです。

実際にベイズ統計解析を行う

§1.準備

実際にRを用いてベイズ統計解析を行ってみたいと思います。
Rでベイズ統計を行うにあたって、環境整備が必要となります。
まず、必要となるものがBUGS/Stanです。
これは、ベイズ推定におけるモデリングを実現するためにどうしても必要となります。
※ モデリングについては前項を参照のこと

今回は導入が比較的容易であり、計算処理が速いと言われているStanを活用しました。

Rstanを利用できる環境にするには、

をインストールすることが必要となります。

§2.そもそもRstanとは何か?

RでRstanを実行してベイズ統計解析を行う前に、ざっくりとRstanに関する説明をすると、

  • stan ・・・ マルコフ連鎖モンテカルロ法を用いたモデリングツール
  • Rstan ・・・ Rを利用してstanへデータ連携するもの

となります。

ここで突然、マルコフ連鎖モンテカルロ法というワードが出てきましたが、大まかに言うと以下のようになります。

  • マルコフ連鎖 ・・・ 前の状態から次の状態が決まる確率過程
    ※例えば、今日の株価は、10年前の株価ではなく、おとといの株価でもなく、昨日の株価から決まるということ。
  • モンテカルロ法 ・・・ ランダムな乱数を使ってシミュレーションを行い確率を算出する

参考までに、マルコフ連鎖モンテカルロ法における乱数発生のアルゴリズムについては以下のものがあります。

  • Metropolis法
  • ギブス・サンプリング法

直接サンプリングが難しい確率分布を持っている場合、サンプリングを近似により行う手法です。
今回はMetropolis法を用いてシミュレーションを行っています。

※Metropolis法のアルゴリズム

metropolis
前章で、『ベイズ統計には主観が伴う』ということに触れました。
これは言い換えると、

パラメータの分布を推定する
※パラメータが確率的に分布するという考えです。

となり、これがリサンプルを作り出すことで真の値に近づいていくであろうという考え方といえます。
parameter_and_data

以下のStanコードを参照してください。

data {
	int<lower=0> N;    // サンプルサイズ
	int<lower=0> y[N]; // 来客者5人当たりの男性客の数
}
parameters {
	real beta; // 全来客数の共通のロジスティック偏回帰係数
	real r[N]; // 男女差
	real<lower=0> s; // 来客者の男女のばらつき
}
model {
	for (i in 1:N)
		y[i] ~ binomial(5, inv_logit(beta+r[i])); // 男性来客率を二項分布で推定
	for (i in 1:N)
		r[i] ~ normal(0,s); // 男女差の事前分布を正規分布で推定
	beta ~ normal(0, 0.001);  // betaを正規分布で推定
	s~uniform(0, 1000);    // sを一様分布で推定
}

これはstan関数と言われるものです。
上から見てみましょう。

dataブロック

π(θ|y)の、θの部分がparametersブロック、yの部分がdataブロックとなります。

data {
	int<lower=0> N;    // サンプルサイズ
	int<lower=0> y[N]; // 来客者5人当たりの男性客の数
}</pre>

サンプルサイズと、目的変数である来客者5人当たりの男性客の数を宣言しています。
与えられた値と、観測値がここに定義されることになります。

parametersブロック

</pre>

parameters {
	real beta; // 全来客数の共通のロジスティック偏回帰係数
	real r[N]; // 男女差
	real<lower=0> s; // 来客者の男女のばらつき
}

推定されるパラメーターを宣言します。
集計における男女差のバラつきが発生するため、推定パラメータに入れています。

modelブロック

</pre>

model {
	for (i in 1:N)
		y[i] ~ binomial(5, inv_logit(beta+r[i])); // 男性来客率を二項分布で推定
	for (i in 1:N)
		r[i] ~ normal(0,s); // 男女差の事前分布を正規分布で推定
	beta ~ normal(0, 0.001);  // betaを正規分布で推定
	s~uniform(0, 1000);    // sを一様分布で推定
}

確率分布での推定を宣言しています。
parameterブロックとmodelブロックにおける主観で事後分布と事前分布を決めています。
男性客の来客は二項分布に従う、事前分布は正規分布に従う、、、などを予め推定しています。
なぜそうなるのか?というと、『従うらしい』です。

§3.Rstanを実行して結果をプロットしてみる

先ほどのstan関数をテキストファイルにコピーし、d.stanとして保存して、
Rコンソールを立ち上げて以下のコマンドを実行します。
なお、今回はレストランの来客者を5人刻みで男女の集計を取ったデータを利用しました。
(ファイル名をcustomer_data.csv としました。)

> require(“rstan”)
> d <- read.csv(“customer_data.csv”)
> head(d)
<< 実行結果 >>
d.csv(“customer_data.csv”)
__id y
1  1 5
2  2 3
3  3 3
4  4 2
5  5 4
6  6 3

続けて以下のコマンドを実行します。

> dat<-list(N=nrow(d),y=d$y)
> d.fit<-stan(file=’d.stanのパス’,data=dat,iter=1000,chains=4)
> plot(d.fit)

fileの設定箇所は、適宜変えてください。
今回は、1000回のシミュレーションを4回セットで実施しました。

■ 以下、実行結果
rstan_result_console

■ 以下、実行結果のプロット図

rstan_result_plot

rを見ると、分布にバラつきがあることは多少はわかると思います。
ただ、もう少し客観的にわかる情報が欲しいですね。
そこで、codaを使って事後の分布関数として見てみましょう。

Rで以下のコマンドを実行します。

> require(“coda”)
> d.fit.coda<-mcmc.list(lapply(1:ncol(d.fit),function(x) mcmc(as.array(d.fit)[,x,])))
> plot(d.fit.coda)

※codaがパッケージに入っていない場合はダウンロードが必要です。

■ 以下、得られる事後分布

posterior_distribution

これが近似的な事後の確率密度分布となります。
(『パラメータがこの辺に存在する確率』を表しています)
これから先に集計していくとsとrの密度より、多少のバラつきはあるものの、男性の割合の方が比較的少なくなりそうです。

男女差の中央値が0付近(若干マイナス)、バラつきの中央値が1.3~1.4として推定されています。

ベイズ推定は、『有限個のデータ』を推定しているため、推定におけるふらつき(曖昧さ)がありますが、
この不確定性が事後分布により表現されています。

まとめ

今回はベイズの推定方法について触れてみました。
パラメーターの推定方法など戸惑う点も多々あるかと思います。
しかし、推定における不確実さもベイズでは考慮されている点はかなり使いやすいと感じました。

今回やっていて気づいた点としては、stanコードでのモデルで収束させるようにするのが思いのほか大変だったことです。
(筆者自身も今回初めてstanに触れたため、なかなか馴染めませんでした)

また、厳密な定義や数式はできる限り省略したので(むしろ省略しすぎましたね)、興味のある方や理解をさらに深めたい方は書籍やWebなどを見て参考してみてください。
以上で今回の記事は以上となります。ご一読ありがとうございました。