需要予測などで、予測の値ではなく分布であると嬉しい場合があります。広告・販促効果分析などで、パラメータ（切片や係数など）が特定の数値ではなく分布であると嬉しい場合があります。
従来の推定方法と同じようなことは、基本ベイズ推定で実施できます。ビジネス活用上、大きな違いの1つは点（1つの値）ではなく分布を手にすることができることです。分布から、平均値や最頻値などを予測値などとすることができます。
ベイズ推定の面白いところは、今手元にあるデータで取り急ぎ計算し、新たなデータを手にしたときに結果を更新していく、というアプローチを取ります。何を更新していくかと言うと、分布を更新していきます。
要は、更新前の分布を、新しいデータを使い、分布を更新していく、ということです。更新前の分布を「事前分布」、更新後の分布を「事後分布」と言います。
そこで問題になってくるのが、最初の事前分布です。
データで更新する前の最初の事前分布なので、何も情報がないということで「無情報分布」という事前分布を利用するケースが多いです。
例えば、一様分布などです。よく利用されます。ただ、推定するパラメータの種類によって無情報分布が変えるのが一般的です。正負の両方を取りうるパラメータのとき正規分布を利用したり、分散などの正の値しかとらない場合はガンマ分布を使用したりします。データ量が増えると、最初の事前分布の影響が小さくなります。
事前分布とデータさえあれば、簡単に事後分布を求めることができるわけではありません。事後分布がいつも簡単に計算できるわけではないからです。
そこで、モンテカルロシミュレーションという乱数を活用した方法が利用されるようになりました。最近では、マルコフ連鎖を利用したモンテカルロシミュレーションであるMCMC（マルコフ連鎖モンテカルロシミュレーション）法を使うのが一般的です。
何となく面倒な気がしますが、ベイズ推定を実現するのがPyMC3というPythonのパッケージです。
今回は、「コイン投げの例を使ってベイズ推定を何となく理解しよう」というお話しをします。今回は、PyMC3は登場しません。
コイン投げ
コイン投げの例で説明を進めていきます。
コインには表と裏があります。コインを投げると、表になったり裏になったりします。それを数字で記録していきます。
- 表：1
- 裏：0
そうすると、次のような0-1のデータが手に入ります。
0010010111011001010101011111000…00101110
100回コイン投げをしたとします。このとき……
- 表；56回
- 裏：44回
……でたとします。
このコインの表の出る確率は、どうでしょうか？
このデータがなければ、0.5と考えるかもしれません。このデータから、0.56（=56/100）と考えるかもしれません。
ベータ分布で考えてみよう
ベータ分布のパラメータを以下のように設定します。
- α = 表の回数56 + 1
- β = 裏の回数44 + 1
ベータ分布に関しては、以下の記事を参考にしてください。
Pythonで描く二項分布とベータ分布
では、ベータ分布Beta(\alpha=57,\beta=45)をPythonで描いてみます。
以下、コードです。
# ベータ分布
## ライブラリー読み込み
from scipy.stats import beta
import numpy as np
import matplotlib.pyplot as plt
## グラフの表示設定
plt.style.use('ggplot') #グラフスタイル
plt.rcParams['figure.figsize'] = [12, 9] # グラフサイズ
## 計算
a = 56+1
b = 44+1
x = np.linspace(0, 1, 100) #x軸
y = beta.pdf(x, a, b) #y軸
## 平均値、最頻値
print("mean = ",a/(a+b))
print("mode = ",(a-1)/(a+b-2))
## グラフ
plt.plot(x, y)
以下、実行結果です。
このベータ分布の横軸は「表の出る確率」で、縦軸が確率密度です。
最も可能性が高い（確率密度の大きい）「表の出る確率」は、最頻値（mode）です。最頻値（mode）は、0.56です。
ベイズ推定とは？
先ず、ベイズの定理を紹介します。
高校生ぐらいで学んだかと思います。高校ぐらいで学んだものとは、ちょっと表記と言うか名称が異なるかもしれません。
以下です。
\displaystyle p(\theta|data)=\frac{p(data|\theta)\cdot p(\theta)}{p(data)}
- 推定対象：\theta（今回の例では、表の出る確率）
- データ：data
- 事前分布（prior）：p(\theta)
- 尤度（likelhood）：p(data|\theta)
- 規格化定数（normalization）：p(data)
- 事後分布（posterior）：p(\theta|data)
これは、得られたデータ（data）を使って、推定対象（\theta）の分布を更新していく式になっています。
更新前の推定対象（\theta）の分布を、事前分布（prior）と言います。
更新後の推定対象（\theta）の分布を、事後分布（posterior）と言います。
最尤法
推定対象（\theta）を推定するとき、従来の統計学では最尤法アプローチをとります。
最尤法は、尤度（likelhood）が最大にする\thetaを探し、それを推定値とするアプローチです。
コインを1回投げたときの尤度p(data|\theta)は以下のようになります。表が「1」で、裏が「0」です。
- p(1|\theta)=\theta
- p(0|\theta)=1-\theta
コインを5回投げたときの尤度p(data|\theta)はどうなるでしょうか？
例えば、表・裏・表・表・裏であった場合、以下のようになります。
\displaystyle p(10110|\theta)=\theta\cdot(1-\theta)\cdot\theta\cdot\theta\cdot(1-\theta)=θ^3\cdot(1-θ)^2
コインを100回投げ、表が56回、裏が44回出た場合には、以下のようになります。
\displaystyle p(data|\theta)=θ^{56} \cdot (1-θ)^{44}
表の出る確率\thetaの推定値は、上の式が最大になる\thetaです。
細かい説明は省きますが、先程登場したベータ分布の最頻値（mode）の0.56す。先程登場したベータ分布は、「表の出る確率」の確率密度です。
ベイズアプローチ
先程示したベイズの定理の数式に、p(data|\theta)=θ^{56} \cdot (1-θ)^{44}を代入してみます。
\displaystyle p(\theta|data)=\frac{p(data|\theta)\cdot p(\theta)}{p(data)}=\frac{θ^{56}\cdot(1-θ)^{44}\cdot p(\theta)}{p(data)}
事後分布p(\theta|data)を得るためには、事前分布p(\theta)と規格化定数p(data)が必要です。
推定対象である\thetaは、今回の例では「表の出る確率」であるため、その事前分布としてベータ分布を想定するのは自然でしょう。繰り返しになりますが、先程登場したベータ分布は「表の出る確率」の確率密度だからです。
ベータ分布は以下です。
\displaystyle f(\theta|\alpha,\beta)=\frac{1}{B(\alpha,\beta)}\theta^{\alpha-1}(1-\theta)^{\beta-1}\\B(\alpha,\beta)=\int_{0}^{1}\theta^{\alpha-1}(1-\theta)^{\beta-1}d\theta
1度もコイン投げをしていない場合、無情報なので、例えば以下のパラメータを設定したベータ分布を事前分布とします。
- α = 表の回数0+1 = 1
- β = 裏の回数0+1 = 1
この場合の事前分布p(\theta)をPythonで描いてみます。
以下、コードです。
## 計算
a = 1
b = 1
x = np.linspace(0, 1, 100) #x軸
y = beta.pdf(x, a, b) #y軸
## 平均値、最頻値
print("mean = ",a/(a+b))
## グラフ
plt.plot(x, y)
以下、実行結果です。
統計学で言うところの、一様分布になります。
この場合の事後分布p(\theta|data)はどうなるのでしょうか？
今回の例では、事前分布p(\theta)をベータ分布と仮定すると、事後分布p(\theta|data)もベータ分布になります。
事前分布がBeta(\alpha,\beta)のもとで、n回施行し表がx回でた場合、事後分布は以下のようになります。
\displaystyle f(\theta|n,x,\alpha,\beta)=\frac{1}{B(x+\alpha,n-x+\beta)}\theta^{x+\alpha-1}(1-\theta)^{n-x+\beta-1}
コイン投げのようなケースは、非常に幸運で簡単に事後分布を導き出すことができます。
通常は、このように簡単にはいきません。p(data)の計算が簡単にいかないからです。
p(data)の計算はベイズ推定アプローチ中で最も難しい部分です。このためのプログラムやライブラリがいくつかあるほどで、そのうちの一つがPyMC3です。
では、この幸運な状況を活用して、無情報下（Beta(\alpha=1,\beta=1)）でコインを100回投げ、表が56回、裏が44回出た場合の事後分布を見てみます。
以下、コードです。
## 計算
a = 56+1
b = 44+1
x = np.linspace(0, 1, 100) #x軸
y = beta.pdf(x, a, b) #y軸
## 平均値、最頻値
print("mean = ",a/(a+b))
print("mode = ",(a-1)/(a+b-2))
## グラフ
plt.plot(x, y)
以下、実行結果です。
これは、Beta(\alpha=57,\beta=45)のベータ分布です。最頻値（mode）は0.56です。
さらにコインを100回投げてみたら、表が46回、裏が54回でした。
この場合の事後分布を見てみます。
事前分布Beta(\alpha=57,\beta=45)の情報下での事後分布を見てみます。
以下、コードです。
## 計算
a = 46+56+1
b = 54+44+1
x = np.linspace(0, 1, 100) #x軸
y = beta.pdf(x, a, b) #y軸
## 平均値、最頻値
print("mean = ",a/(a+b))
print("mode = ",(a-1)/(a+b-2))
## グラフ
plt.plot(x, y)
以下、実行結果です。
まとめ
今回は、「コイン投げの例を使ってベイズ推定を何となく理解しよう」というお話しをしました。
今回は、PyMC3は登場しませんでした。コイン投げの場合、登場するまでもない、という感じでしたが、通常はそうは行きません。
次回は、「コイン投げの例を使ってPyMC3ベイズ推定を何となく理解しよう」というお話しをします。
