Juliaでコイントスのパラメータの事後分布を推定する

Juliaで作って学ぶベイズ統計学という本を読んでいます。

用語の確認

• 同時確率
  • 複数の事象が同時に起きる確率
  • $P(X, Y)$
• 周辺確率
  • 同時確率から複数の事象(確率変数)を除いた確率
  • $P(X) = \sum_Y P(X, Y)$
• 条件付き確率
  • ある事象が起こる確率における、ある事象の確率(言語化むずい)
  • Yが起こる事象のうち、Xが起こる事象の確率
  • $P(X|Y) = \frac{P(X,Y)}{P(Y)}$
• 事前分布
  • パラメータを生成する確率分布
• 事後分布
  • 観測結果という制約が入った前提で、パラメータを生成する尤もらしい確率分布

生半可な考えだが、事前分布は時間に順行しているが、事後分布は時間に逆行した視点だと解釈している。事後分布は、観測結果という強い証拠を持ち合わせている。

本題

Aさんは、コイントスで表が出る確率をランダムに決める。その確率をもとに、 $N$ 回コイントスを行う。Bさんは、コイントスの結果から表が出る確率を当てたい。

これをかしこまった言い方にして、解いていく。

パラメータ $\mu$ は一様分布からサンプリングされる。

ベルヌーイ分布 $\text{Ber}(\mu)$ から $N$ 回試行された観測結果 $\bm{X} = \{ x_1, x_2, \dots , x_N\}$ から、パラメータ $\mu$ の事後分布を推定する。

$$\mu \sim \text{Uniform}(\mu | 0, 1) \\ x_n \sim \text{Ber}(x_n | \mu)$$

観測結果が起こるような同時分布は次のように書ける。

$$\begin{align} p(x_1, x_2, \dots , x_N, \mu) &= p(\mu) \prod^N_{n=1} p(x_n | \mu)\\ &= p(\mu) p(\bm{X} |\mu) \end{align}$$

求めたい観測結果 $\bm{X}$ が与えられたときのパラメータ $\mu$ の事後分布は、条件付き分布の定義から次のようになる。

$$p(\mu | \bm{X}) = \frac{p( \bm{X} , \mu ) }{p(\bm{X} )}$$

分子については、ベルヌーイ分布と一様分布の同時確率として与えられるが、分母については計算が必要になる。

分母 $p(\bm{X})$ は、分子をパラメータ $\mu$ について周辺化して計算できる。$p(\bm{X})$ を周辺尤度と呼ぶ。

$$p(\bm{X}) = \int p( \bm{X} , \mu ) d\mu$$

まとめると、求めたい観測結果 $\bm{X}$ が与えられたときのパラメータ $\mu$ の事後分布は次のように求まる。

$$p(\mu | \bm{X}) = \frac{p( \bm{X} , \mu ) }{ \int p( \bm{X} , \mu ) d\mu}$$

これをJuliaで数値積分をすることによって具体的な計算を行う。

# 観測結果
X_obs = [0, 0, 0, 1, 1, 1, 1, 1]

# 同時確率分布 p(X, μ), 分子に相当
p_joint(X, μ) = prod(pdf.(Bernoulli(μ), X)) * pdf(Uniform(0, 1), μ)

# 周辺尤度を計算する関数を返す
function approx_integration(μ_range, p)
    Δ = μ_range[2] - μ_range[1]
    X -> sum([p(X, μ) * Δ for μ in μ_range])
end


# ベルヌーイ分布なのでμの積分範囲は0 ~ 1
μ_range = 0:1e-2:1

# 周辺尤度, 分母に相当
p_marginal = approx_integration(μ_range, p_joint)

# 事後分布
postprior(μ) = p_joint(X_obs, μ) / p_marginal(X_obs)

plot(postprior, label="", xlabel="μ")

このコードを実行すると、次のようなパラメータ $\mu$ の事後分布のグラフが出てきます。

参考文献

• Juliaで作って学ぶベイズ統計学 | 講談社