この記事で扱う内容
ベイズ推定の概要
- 統計モデリングとは
- パラメタの推定手法としてのベイズ推定の利点
ベルヌーイ分布モデルのベイズ的解釈
- ベルヌーイ分布モデルとは
- ガンマ分布を事前分布に用いた場合の事後分布と予測分布
Pythonによる数値実験
- 事後分布が最尤推定解付近に収束していく様子
- 真の分布と予測分布間のKL ダイバージェンス
参考文献

この記事で扱う内容

　本稿は、統計モデリングの基本中の基本であるベルヌーイ分布モデルを題材に、初学者が統計モデリングにおけるベイズ推定の有用性を体感することを目的としています。統計モデリングでは、観測データはある確率分布（以降は"真の確率分布"と呼ぶことにします。）に従い観測される確率変数であると仮定し、観測データの生成過程を確率モデルの組み合わせで表現することを考えます。中でもベルヌーイ分布は、「バク転が成功するか失敗するか」「コインを投げて表が出るか裏が出るか」「ホークスが今日の試合に勝つか負けるか」などの同時には起こり得ない２つの事象を表現する際に使用される確率モデルです。本稿では、ベイズ推論について概要を簡単に説明した後、ベルヌーイ分布モデルのベイズ的解釈について述べます。その後、Pythonによる簡単なデモンストレーションを行い、観測データが増加するにつれてパラメタの事後分布が真の分布に収束していく様子を確認します。

ベイズ推定の概要

　まずは、統計モデリングとベイズ推定の一般論を確認しましょう。全体感重視なざっくり目の説明なので、式変形などの細部を詰めたい人は末記のリンクや標準的な数理統計学の教科書を参照してください。

統計モデリングとは

　統計モデリングとは一言で言えば、確率変数 ${X}$ を生成する真の確率分布 ${q(X)}$ をモデル設計者が恣意的に設計した確率分布 ${p(X)}$ により近似的に表現するための営みのことです。モデル設計者は、真の確率分布 ${q(X)}$ の近似モデル ${p(X)}$ を確率モデルの組み合わせとして構成し、真の確率分布 ${q(X)}$ から得られた観測データ $\mathcal{D}$ を利用してモデルのパラメタ ${\theta}$ に関する推論を行います。ここでパラメタ $\theta$ とはモデルの表式に陽に現れる変量のことであり、パラメタの値がモデル分布の形状を決定します。つまり同じ確率分布のクラスであっても、パラメタの値が異なれば分布の形状は全く異なるのです。例えば、下の図の３種のヒストグラムは全てポアソン分布と呼ばれる確率分布から生成されたものですが、パラメタの値が異なるために分布の形状が異なっています。 f:id:snow67675476:20190910031755p:plain

以上を踏まえれば、モデル設計者はモデルの確率分布のクラスを定めることに加えて、モデル ${p(X)}$ が真の分布 ${q(X)}$ に少しでも近づくようにパラメタ $\theta$ に関して何らかの推論を行う必要があることは想像に難くないでしょう。即ち、統計モデリングを正当に行うためには、「モデル ${p(X, \theta)}$ の確率分布クラスの選定」と「モデルのパラメタ $\theta$ に関する推論」の２つの関門を突破する必要があるということが分かります。 *1

パラメタの推定手法としてのベイズ推定の利点

　前節では、統計モデリングの課題は「モデル ${p(X, \theta)}$ の確率分布クラスの選定」と「モデルのパラメタ $\theta$ に関する推論」の２部に分けられることを見ました。本稿の主要テーマであるベイズ推定は「モデルのパラメタ $\theta$ に関する推論」を行うためのメジャーな手法の１つです。この節では、パラメタ推定におけるベイズ推定の特徴を他の類似手法と比較しながら見ていきます。

　さて、モデル ${p(X, \theta)}$ を表す確率分布のクラスを"適切に"仮定することが出来たとしましょう。次にモデル設計者が乗り越えなくてはならない関門は、モデル ${p(X, \theta)}$ ができる限り真の分布 ${q(X)}$ の特徴を捉えることが出来るようにパラメタ $\theta$ に関して妥当な推論を行うことです。パラメタに関する推論の方針は大別すると３パターンあり、１つ目はデータ既知の下で尤度関数 ${p(D|\theta)}$ を最大化するパラメタを点推定する最尤推定、２つ目はデータ既知の下でパラメタの事後分布 ${p(\theta|D)}$ を最大化するパラメタを点推定するMAP推定、そして３つ目は尤度関数 ${p(D|\theta)}$ のパラメタの事後分布 ${p(\theta|D)}$ に関する期待値を予測分布 ${p(X|\mathcal{D})}$ とするベイズ推定です。以降では、各手法の特徴を簡単に見ていきます。

　最尤法とMAP推定は目的関数を最大化するという意味で"最適な"パラメタ ${\hat{\theta}}$ が一意に定める点推定法であり、モデル ${p(X|\hat{\theta})}$ を真の確率分布 ${q(X)}$ の近似分布とします。目的関数の最大化という基準でパラメタを点推定する方針は一見合理的に思われますが、ここで言う"最適"なパラメタとはあくまで「観測データに対して最も適合する」という意味での"最適解"に過ぎず、未知のデータに対する汎化性能を保証するものではありません。実際、最尤推定やMAP推定はモデルのパラメタ数や観測データの次元数が増大するにつれて過学習を起こしてしまうリスクが著しく大きくなることが知られており、実用上はパラメタを点推定することに対して慎重になるべきです。一方で、３つ目のベイズ推定は"最適な"パラメタを一意に定めることはせずに、事後分布 ${p(\theta|D)}$ や予測分布 ${p(X|\mathcal{D})}$ といった確率分分布の形式で推論結果を得ます。従って、ベイズ推定ではパラメタに関する推論の"不確実性"を定量的に表現することが出来るという利点があります。また、予測分布 ${p(X|\mathcal{D})}$ の計算にはパラメタの定義域上の確率分布である事後分布 ${p(\theta|D)}$ の情報を最大限に活用するため、"最適解"以外の情報を捨て去ってしまう最尤法やMAP推定と比較して、過学習を起こしにくいという利点もあります。予測分布の具体的な計算手順については、次項の具体例を見ながら確認していきます。

datachemeng.comwww.hellocybernetics.tech

ベルヌーイ分布モデルのベイズ的解釈

　まずは、ベルヌーイ分布モデルの概要を確認します。その後、ベルヌーイ分布モデルのベイズ的解釈において、事前分布にベータ分布を用いればパラメタの事後分布と予測分布が解析的に計算できることを確認します。

ベルヌーイ分布モデルとは

　例えば、確率 ${\mu}$ でバク転を成功できる人が ${N}$ 回バク転を実施した時、「バク転が成功した時に ${1}$ 、失敗した時に ${0}$ 」となる確率変数 ${X}$ の実現値 ${x \in \{0, 1\}}$ の系列を観測データ ${{\cal D}=\{x_1, x_2, \cdots , x_N\}\subset \{0, 1\}^ N}$ として記録していくことを考えます。この観測データ ${\cal D}$ を用いて、バク転の成功確率 $\mu$ に関して何らかの知見を得たいとします。ここでモデル設計者が、 $\forall{i} \in \{ 1, 2, \cdots, N \}$ に対して「 ${i}$ 回目の試行に対応する確率変数 ${X_i}$ は確率 $\mu$ で ${1}$ 、確率 ${1-\mu}$ で ${0}$ となる」*2との仮説を立てたとしましょう。すると確率変数 ${X_i}$ の確率分布関数は以下のように表現することが出来ます。

$\begin{align} p(X_i|\mu) &= \mu^ {X_i}(1-\mu)^ {1-X_i} \Leftrightarrow X_i \sim \mathrm{Bernoulli}(\mu) \tag{1}\\ \end{align}$

このように、離散２値変数 ${x \in \{0, 1\}}$ 上の確率分布であり、分布のパラメタ ${\mu}$ の値に応じて ${0}$ と ${1}$ の生成確率が定まるような分布をベルヌーイ分布と呼びます。上記の例のように、ベルヌーイ分布を用いて現象をモデリングすることをベルヌーイ分布モデルと呼ぶことにします。

　さて、バク転の成功確率 $\mu$ はベルヌーイ分布モデルのパラメタなので、観測データ ${\cal D}$ を用いて１）最尤推定やMAP推定によりパラメタを点推定する２）ベイズ推定によりパラメタの事後分布と予測分布を計算するなどの方針で、バク転の成功確率 $\mu$ に関して知見を得る事が出来ると期待できます。今回はベイズモデリングに焦点を当てるため、後者のアプローチによるパラメタ推定を考えたいと思います。

ガンマ分布を事前分布に用いた場合の事後分布と予測分布

　一般に、ベイズモデリングを行う際は、真のデータ生成分布 ${q(X)}$ の近似表現として個々の観測データの生成モデル ${p (X, \theta)}$ を仮定し、全体の観測データ ${\cal D}$ は真のデータ生成分布 ${q(X)}$ からのランダムなサンプリングとして得られていると考えます。ベイズの定理によれば、全体の観測データの生成モデル ${p (\cal{D}, \theta)}$ は 尤度関数 ${p ({\cal D}| \theta)}$ と事前分布 ${p (\theta)}$ の積に分解出来るので、モデル設計者は尤度関数 ${p ({\cal D}| \theta)}$ と事前分布 ${p (\theta)}$ の確率分布のクラスを仮定する必要があります。今回は「バク転が成功するか否か」という2値の離散データのモデリングを考えるので、各データのモデル ${ \left\{p (X_i| \theta)\right\}_{i=1}^{N} }$ としては全て同一のベルヌーイ分布 ${\mathrm{Bernoulli}(\mu)}$ と仮定するのが一般的でしょう。それでは、パラメタの事前分布 ${p (\theta)}$ の確率分布クラスについてはどのような基準に基づいて設定すれば良いのでしょうか？結論から言うと、パラメタの事前分布 ${p (\theta)}$ としては、 ${\mu \in (0, 1)}$ なる実数を生成する分布であるベータ分布 ${\mathrm{Beta}(a, b)}$ *3を選択することが多いです。これはベータ分布がベルヌーイ分布の共役事前分布であることに由来します。共役事前分布とは事後分布と事前分布が同じ確率分布クラスになるような事前分布のことです。事前分布と事後分布の確率分布が一致していると、事後分布や予測分布を導出するための計算が煩雑になりにくいという利点があるため、特別な理由がない限りは事前分布は尤度関数の共役事前分布を選択することが多いようです。以降では例によって、尤度関数 ${p ({\cal D}| \theta)}$ としてベルヌーイ分布、事前分布 ${p (\theta)}$ としてベータ分布を仮定して議論を進めることにします。事前分布の選択に関しては少々乱暴な議論に思えるかもしれませんが、観測データが増加するにつれて事前分布が予測に与える影響は少なくなっていくので、計算効率を優先したモデル設計を行う方が結果としてうまくいくことが多いです。

　少し前置きが長くなりましたが、ここからはベルヌーイ-ベータモデルの事後分布と予測分布を導出していこうと思います。以下では、観測データ ${{\cal D}=\{x_1, x_2, \cdots , x_N\}\subset \{0, 1\}^ N}$ が独立に得られていると仮定します。まずは、パラメタの事後分布に関してベイズの定理と積の法則を用いて式変形を行います。

$\begin{align} p(\theta|{\cal D}) &= \frac {\displaystyle p({\cal D}|\theta)p(\theta)} {\displaystyle p({\cal D})} \\ &= \frac {\displaystyle \left\{\prod_{i=1}^ {N} p ( x_i|\theta)\right\}p(\theta)} {\displaystyle p({\cal D})} \\ &\propto \left\{ \prod_{i=1}^ {N} p ( x_i|\theta)\right\}p(\theta) \tag{2} \\ \end{align}$

　パラメータの事後分布 ${p ( \theta|{\cal D})}$ に対して、ベイズの定理と積の法則を適用することで（２）式を得ます。なお、分布の形状に影響を与えるのはパラメタ ${\theta}$ に依存する部分のみなので、２行目から３行目で正規化係数 ${p({\cal D})}$ を省略しています。ここまでの変形は観測データ ${\cal D}$ の独立性のみを仮定しており、モデル ${p (X, \theta)}$ の設計によらず成立します。以降は簡単のために（２）式に対数をとって計算を進めていくことにします。

$\begin{align} \ln {p(\theta|{\cal D})} &= \sum_{i=1}^ {N} \ln {p(x_i|\theta)}+\ln {p(\theta)} + \mathrm{const.} \\ &= \sum_{i=1}^ {N} X_i \ln {\mu}+\sum_{i=1}^ {N} (1-x_i) \ln {(1-\mu)} +(a-1) \ln {\mu}+(b-1)\ln {(1-\mu)} + \mathrm{const.} \\ &= (\sum_{i=1}^ {N} x_i +a-1)\ln {\mu}+(N-\sum_{i=1}^ {N} X_i +b-1)\ln{(1-\mu)} + \mathrm{const.} \tag{3} \\ \end{align}$

${p(x_i|\theta)=\mathrm{Bernoulli}(\mu)=\mu^ {x_i}(1-\mu)^ {1-x_i}}$ と ${p(\mu)=\mathrm{Beta}(a, b)=C_{B}(a, b)\mu^ {a-1}(1-\mu)^ {b-1} }$ を代入して変形すれば（３）式を得ます。ハイパーパラメタ $a, b \in \mathbb{R}^ {+}$ は、モデル設計者が定義域内で適切に定めた"定数"*4であり、規格化因子 $C_B(a, b)$ *5もまたハイパーパラメタの関数として定まる"定数"です。 ${\mathrm{Beta}(a, b)}$ の対数を取ったものを計算してみると以下のようになります。

$\begin{align} \ln {\mathrm{Beta}(a, b)} &= (a-1) \ln{\mu}+(b-1)\ln{(1-\mu)}+\mathrm{const.} \tag{4} \\ \end{align}$

（３）式と（４）式を比較すると、

$\begin{align} p(\mu|{\cal D})&=\mathrm{Beta}(\hat {a}, \hat {b})\tag{5} \\ \end{align}$

$\begin{align} \left\{ \begin{array}{l} \hat{a} = \sum_{i=1}^ {N} x_i + a \tag{6}\\ \hat{b} = N-\sum_{i=1}^ {N} x_i+b \\ \end{array} \right. \end{align}$

となることが分かり、パラメタの事後分布が事前分布と同じくベータ分布になることを解析的に導出することが出来ました。ただし、新しいベータ分布のパラメタ $\hat{a}$ と $\hat{b}$ は観測データ ${\cal D}$ に依存する形でベイズ更新されているので、分布の形状は事前分布とは異なることに注意してください。*6 次に、予測分布 ${p(X_\mathrm{new}|\cal{D})}$ を導出していきます。観測データ ${\cal D}$ が所与である場合の予測分布 ${p(X_\mathrm{new}|\cal{D})}$ は、*7 ベイズ更新後のモデル ${p (X_{\mathrm{new}}, \mu|{\cal D})}$ の確率変数 ${X_\mathrm{new}}$ に関する周辺分布なので、モデルのパラメタ $\mu$ を積分除去することにより導出することが出来ます。

$\begin{align} p ( X_\mathrm{new}|{\cal D}) &=\int_{ } p (X_{\mathrm{new}}, \mu|{\cal D}) d\mu \\ &=\int_{ } p (X_{\mathrm{new}}|\mu)p(\mu|{\cal D}) d\mu \\ &=\int_{ } \mathrm{Bernoulli}(\mu)\mathrm{Beta}(\hat{a}, \hat{b}) d\mu \\ &=\int_{ } \mu^ {X_\mathrm{new}}(1-\mu)^ {1-X_\mathrm{new}}C_{B}(\hat{a}, \hat{b})\mu^ {\hat{a}-1}(1-\mu)^ {\hat{b}-1} d\mu \\ &=C_{B}(\hat{a}, \hat{b})\int_{ } \mu^ {X_\mathrm{new}+\hat{a}-1}(1-\mu)^ {1-X_\mathrm{new}+\hat{b}-1} d\mu \tag{7}\\ \end{align}$

１行目から２行目はベイズの定理を用いて変形しています。２行目の式は、まさに尤度関数の事後分布に関する期待値の定義そのものになっていますね。ここで（７）式の積分の部分はベータ分布の定義式を当てはめれば、

$\begin{align} \int_{ } \mu^ {X_\mathrm{new}+\hat{a}-1}(1-\mu)^ {1-X_\mathrm{new}+\hat{b}-1} d\mu= \frac {\displaystyle 1}{\displaystyle C_B(X_\mathrm{new}+\hat{a}, 1-X_\mathrm{new}+\hat{b})} \tag{8} \\ \end{align}$

となるので、予測分布は

$\begin{align} p ( X_\mathrm{new}|{\cal D}) &= \frac {\displaystyle C_B(\hat{a}, \hat{b})}{\displaystyle C_B(X_\mathrm{new}+\hat{a}, 1-X_\mathrm{new}+\hat{b})} \tag{9} \\ &=\frac {\displaystyle \Gamma(\hat{a}+\hat{b})\Gamma(X_\mathrm{new}+\hat{a})\Gamma(1-X_\mathrm{new}+\hat{b})}{\displaystyle \Gamma(\hat{a})\Gamma(\hat{b})\Gamma(\hat{a}+\hat{b}+1)} \tag{9} \\ \end{align}$

と書くことができます。ここでは、規格化因子の具体的な表示が、

$\begin{align} C_B(a, b)&= \frac {\displaystyle \Gamma(a+b)}{\displaystyle \Gamma(a)\Gamma(b)} \tag{10} \\ \end{align}$

となることを用いました。予測分布の標識はかなり複雑に見えますが、ガンマ関数 $\Gamma(\cdot)$ が満たす等式

$\begin{align} {\Gamma(x+1)=x\Gamma(x)} \tag{11} \\ \end{align}$

を用いて分解すれば、 ${p(X_\mathrm{new}=1|\mathcal{D})}$ は次のように簡潔な形で書くことができます。

$\begin{align} p ( X_\mathrm{new}=1|{\cal D}) &= \frac {\displaystyle \Gamma(\hat{a}+\hat{b})\Gamma(1+\hat{a})\Gamma(\hat{b})}{\displaystyle \Gamma(\hat{a})\Gamma(\hat{b})\Gamma(\hat{a}+\hat{b}+1)} \\ &=\frac {\hat{a}}{\hat{a}+\hat{b}} \tag{12} \\ \end{align}$

${p(X_\mathrm{new}=1|\mathcal{D})}$ についても同様にして、

$\begin{align} p ( X_\mathrm{new}=0|{\cal D}) &= \frac {\displaystyle \Gamma(\hat{a}+\hat{b})\Gamma(\hat{a})\Gamma(1+\hat{b})}{\displaystyle \Gamma(\hat{a})\Gamma(\hat{b})\Gamma(\hat{a}+\hat{b}+1)} \\ &=\frac {\hat{b}}{\hat{a}+\hat{b}} \tag{13} \\ \end{align}$

と簡潔に書くことができます。（１２）式、（１３）式を１つの式にまとめると、

$\begin{align} p ( X_\mathrm{new}|{\cal D}) &= \left(\frac {\displaystyle \hat{a}}{\displaystyle \hat{a}+\hat{b}}\right)^{X_\mathrm{new}}\left(\frac {\displaystyle \hat{b}}{\displaystyle \hat{a}+\hat{b}}\right)^{1-X_\mathrm{new}} \\ &=\left(\frac {\hat{a}}{\hat{a}+\hat{b}}\right)^{X_\mathrm{new}}\left(1-\frac {\hat{a}}{\hat{a}+\hat{b}}\right)^{1-X_\mathrm{new}} \\ &=\mathrm{Bernoulli}\left(\frac {\displaystyle \hat{a}}{\displaystyle \hat{a}+\hat{b}} \right) \tag{14} \\ \end{align}$

となり、予測分布はベルヌーイ分布となることを解析的に導くことができました。予測分布のパラメタは ${\hat{a}}$ と ${\hat{b}}$ の有理式になっていますが、各々は（６）式の形で書けるのでした。ここで、 ${ \sum_{i=1}^{N}x_i}$ はバク転の成功回数を意味することに注意すれば、バク転の成功回数が多ければ多いほど、 ${\hat{a}}$ の値もそれに比例して大きくなることが分かります。一方で、 ${N-\sum_{i=1}^{N}x_i}$ はバク転の失敗回数を意味するので、バク転の失敗回数が多ければ多いほど、 ${\hat{b}}$ の値もそれに比例して大きくなります。従って予測分布のパラメタは、事前分布のハイパーパラメタである ${a}$ と ${b}$ をオフセットとする観測データ $\mathcal{D}$ 中のバク転の成功確率を表していると解釈できます。これは、「被験者がバク転を成功すればするほど、将来もバク転を成功する確率が高いであろう」と期待することに他ならず、直感的にも推論の妥当性が理解できるでしょう。

　この節では、ベルヌーイ-ベータモデルのパラメタの事後分布と予測分布を解析的に導出し、その直感的な意味について考察しました。一般には、モデルが複雑になればなるほどパラメタの事後分布と予測分布を解析的に求めることが困難になります。そのような場合には、事後分布や予測分布の導出はMCMC(Markov Chain Monte Carlo)や変分推論などの近似手法に頼ることになります。解析解が求まらない場合の近似手法については、以下のリンクを参照して下さい。

machine-learning.hatenablog.com

Pythonによる数値実験

最後に、Pythonによる数値実験を行います。以下のソースコードはAnaconda+TensorFlow Probabilityの環境で動作します。残念ながらTensorFlow ProbabilityはPython3.7系以降にまだ対応していないので、環境構築が面倒な人はGoogle Colaboratoryを使うことをお勧めします。

事後分布が最尤推定解付近に収束していく様子

下図では、観測データ数が増加していくにつれて事後分布の分散が小さくなっていく様子をプロットしています。データ数０の場合は事前分布を表しています。データ数が増えるにつれて、徐々に分布の中心が最尤推定解付近に収束すると同時に分散も小さくなっていくことが分かります。ソースコードも示しますので、ハイパーパラメタやバク転の試行回数の値を変更して、色々と試してみるとより理解が深まると思います。

f:id:snow67675476:20190916193015p:plain


%matplotlib inline
from IPython.core.pylabtools import figsize
import numpy as np
from matplotlib import pyplot as plt
figsize(15, 12)

import scipy.stats as stats

"""事前分布"""
dist = stats.beta
a = 1.0
b = 10.0

"""バク転の試行回数"""
n_trials = [0,10, 20, 30, 100, 200, 300, 500, 700, 1000]

"""バク転の成功確率は50パーセントとし、
ベルヌーイ分布からサンプリングして観測データを作成する
これは真の分布なので本来は未知である"""
mu = 0.5
data = stats.bernoulli.rvs(mu, size=n_trials[-1])

x = np.linspace(0, 1, 100)

for k, N in enumerate(n_trials):
    sx = plt.subplot(len(n_trials)//2, 2, k+1)
    plt.xlabel("$p$, probability of success") if k in [0, len(n_trials)-1] else None
    plt.setp(sx.get_yticklabels(), visible=False)
    
    """N回目までのバク転の成功回数"""
    success = data[:N].sum()
    
    """パラメタの事後分布"""
    posterior = dist.pdf(x, success+a, N-success+b)
    
    plt.plot(x, posterior, label= '%d trial, \n %d success'%(N, success))
    plt.fill_between(x, 0, posterior, color='#348ABD', alpha=0.4)
    plt.vlines(mu, 0, 4, color='k', linestyles='--', lw=1)
  
    leg = plt.legend()
    leg.get_frame().set_alpha(0.4)
    plt.autoscale(tight=True)
  
plt.suptitle('Baysian updating of posterior probabilities', y=1.02, fontsize=14)
plt.tight_layout()

真の分布と予測分布間のKL ダイバージェンス

次に、真の分布 ${q(X)}$ と予測分布 $p(X|\mathcal{D})$ 間のKL ダイバージェンスを観測データ $\mathcal{D}$ に対してプロットしてみます。 KLダイバージェンスは２つの確率分布間の"距離"を測る尺度であり、値が0に近いほど２つの確率分布が良く似ていると言えます。KLダイバージェンスの概要については以下のリンクを参考にしてください。

qiita.com

下図を見ると観測データ数が増加するにつれたKLダイバージェンスが０に漸近していく様子が分かりますね。今回は非常に単純なモデルということもあり、観測データ数が100を超えたあたりからは、予測分布が真の分布をほぼ表現できているといえます。

f:id:snow67675476:20190916192456p:plain


 from math import log

"""真の分布のエントロピー"""
entropy = -mu*log(mu)-(1-mu)*log(1-mu)

KL = []
for N in range(n_trials[-1]):
    """N回目までのバク転の成功回数"""
    success = data[:N].sum()
    
    """"予測分布のパラメタ"""
    mu_hat = (success+a)/(N+a+b)
    
    """KL ダイバージェンス"""
    KL.append(-entropy-mu*log(mu_hat)-(1-mu)*log(1-mu_hat))
    
x = np.arange(n_trials[-1]) 
plt.plot(x, KL,  color='#348ABD', alpha=0.9)
plt.xlim(0, 1000)

p = plt.hlines([0], 0, 1000, "red", linestyles='dashed')  
plt.title('KL divergence between of q(X) and p(X|D)',  fontsize=20)
plt.xlabel("# of trials", fontsize=20)
plt.ylabel("KL divergence", fontsize=20)

参考文献

ベイズ推論による機械学習入門　著　須山敦志
Probabilistic-Programming-and-Bayesian-Methods-for-Hackers by Cameron Davidson-Pilon

*1:モデルのパラメタ $\theta$ は観測データ $\mathcal{D}$ を用いて推定すべき対象なので確率変数として扱います。以降はモデルを表記する際は ${p(X, \theta)}$ のようにパラメタ $\theta$ を引数に明記することにします。

*2:確率変数 ${X_i}$ が従う確率分布 ${p (X_i| \mu)}$ のパラメタ ${\mu}$ が ${i}$ の値に依存しないことが、各 ${X_i}$ が同分布からサンプリングされていると仮定することに対応します。

*3: ${a}$ と ${b}$ は正の実数であり、パラメタの分布のパラメタということでハイパーパラメタと呼ばれています。この値はモデル設計者が適切に値を決める必要があります。

*4:今回のモデルにおいては観測データ ${\cal D}$ を用いて推定を行う対象としない、即ち確率変数ではないという意味で"定数"と表現しています。

*5:具体的な式は式(１０)で与えられます。

*6: ${\hat{a}}$ と ${\hat{b}}$ は確率変数 ${X}$ に依存するので確率変数です。

*7:観測データが所与でない場合も、パラメタの事前分布を用いて予測分布を計算することは可能です。

機械学習チョットデキルようになりたい

脱サラ大学院生が機械学習を研究してみた。機械学習や統計学を中心に、日々勉強したことを投稿していきます。

ベルヌーイ分布モデルで体験するベイズ推定