Juliaを使ってみる(10/22)：厳密解法によるベルヌーイモデル

今回やりたいこと

　確率試行のシミュレーションをJuliaで扱ってみる。

Juliaで作って学ぶベイズ統計学 (KS情報科学専門書)

を参照する。

前回

power-of-awareness.com

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))")

次回

power-of-awareness.com

*1:直下のfor文内におけるの実装が参照テキスト(第1刷)P.144では逆になっている。同ページ下部の予測に関する実装は直下のとおりになっているため、本ブログの内容が正しいと考える。