「大人の教養・知識・気付き」を伸ばすブログ

一流の大人(ビジネスマン、政治家、リーダー…)として知っておきたい、教養・社会動向を意外なところから取り上げ学ぶことで“気付く力”を伸ばすブログです。目下、データ分析・語学に力点を置いています。今月(2022年10月)からは多忙につき、日々の投稿数を減らします。

MENU

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

今回やりたいこと

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

を参照する。

1. 統計モデリングと推論

1.1 厳密解法

 ベルヌーイ・ベータモデルは共役事前分布を用いたシンプルなモデルであり、厳密な事後分布を解析的に求めることができることが知られている。ベータ分布Beta(\alpha,\beta)


\begin{aligned}
\alpha&=1+\displaystyle{\sum_{n=1}^{N}x_n},\\
\beta&=1+N-\displaystyle{\sum_{n=1}^{N}x_n}
\end{aligned}


\muの事後分布である。これに基づき事後分布の密度関数をプロットしてみる*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文内における\alpha,\betaの実装が参照テキスト(第1刷)P.144では逆になっている。同ページ下部の予測に関する実装は直下のとおりになっているため、本ブログの内容が正しいと考える。

プライバシーポリシー お問い合わせ