今回やりたいこと
確率試行のシミュレーションをJuliaで扱ってみる。
を参照する。
1. 統計モデリングと推論
1.1 厳密解法
ベルヌーイ・ベータモデルは共役事前分布を用いたシンプルなモデルであり、厳密な事後分布を解析的に求めることができることが知られている。ベータ分布
がの事後分布である。これに基づき事後分布の密度関数をプロットしてみる*1。
using Distributions, PyPlot, LinearAlgebra # グラフの軸やタイトルを設定する関数 function set_options(ax, xlabel, ylabel, title; grid = true, gridy = false, legend = false) ax.set_xlabel(xlabel) ax.set_ylabel(ylabel) ax.set_title(title) if grid if gridy ax.grid(axis = "y") else ax.grid() end end legend && ax.legend() end ### 観測値 X_obs1 = [0,0,0,1,1] X_obs2 = [0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1] fig, axes = subplots(1,2,sharey = true, figsize = (10,4)) mus = range(0, 1, length = 100) for (i, X_obs) in enumerate([X_obs1, X_obs2]) # 厳密な事後分布:ベータ分布 alpha = 1.0 + sum(X_obs) beta = 1.0 + length(X_obs) - sum(X_obs) d = Beta(alpha, beta) # 事後分布を可視化 axes[i].plot(mus, pdf.(d, mus)) set_options(axes[i], "mu", "density", "exact posterior (X_obs#(i))") end

# 予測 function prediction(X_obs) alpha = 1.0 + sum(X_obs) beta = 1.0 + length(X_obs) - sum(X_obs) alpha/(alpha + beta) end println("$(prediction(X_obs1)), $(prediction(X_obs2))")
*1:直下のfor文内におけるの実装が参照テキスト(第1刷)P.144では逆になっている。同ページ下部の予測に関する実装は直下のとおりになっているため、本ブログの内容が正しいと考える。