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

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

MENU

Juliaを使ってみる(12/22):数値積分を用いたシミュレーション

今回やりたいこと

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

を参照する。

1. 数値積分を使った線形モデルの学習

1.1 理論の準備

 条件付き分布の定義から、学習用のデータセット\boldsymbol{X}および\boldsymbol{Y}が与えられた後の母数の事後分布は



\begin{aligned}
P\{\boldsymbol{w}|\boldsymbol{X},\boldsymbol{Y}\}=\displaystyle{\frac{P\{\boldsymbol{Y}|\boldsymbol{w},\boldsymbol{X}\}P\{\boldsymbol{w}\}}{P\{\boldsymbol{Y}|\boldsymbol{X}\}}}
\end{aligned}


で与えることができる。ここで



\begin{aligned}
P\{\boldsymbol{Y}|\boldsymbol{X}\}=\displaystyle{\int P\{\boldsymbol{Y}|\boldsymbol{w},\boldsymbol{X}\}P\{\boldsymbol{w}\}}d\boldsymbol{w}
\end{aligned}


である。したがって上記の積分計算ができるならば、事後分布の密度関数が計算できる。

1.1.1 シミュレーション例
using Distributions
using PyPlot


# 同時分布
p_joint(X, Y, w₁,w₂) = prod(pdf.(Normal.(w₁.*X.+w₂,σ), Y)) * 
                       pdf(Normal(μ₁,σ₁), w₁) * pdf(Normal(μ₂,σ₂),w₂)


# 数値積分
function approx_integration_2D(w_range, p)
    Δ = w_range[2] - w_range[1]
    (X, Y) -> sum([p(X, Y, w₁, w₂) * Δ^2 for w₁ in w_range, w₂ in w_range])
end


### シミュレーション
X_obs = [-2, 1, 5]
Y_obs = [-2.2, -1.0, 1.5]

σ = 1.0
μ₁ = 0.0
μ₂ = 0.0
σ₁ = 10.0
σ₂ = 10.0

# wの積分範囲
w_range = range(-5,5,length=500)


# 数値積分を実行
p_marginal = approx_integration_2D(w_range, p_joint)

p_marginal(X_obs, Y_obs)



using Distributions
using PyPlot

# グラフの諸設定を行う関数
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


# 同時分布
p_joint(X, Y, w₁,w₂) = prod(pdf.(Normal.(w₁.*X.+w₂,σ), Y)) * 
                       pdf(Normal(μ₁,σ₁), w₁) * pdf(Normal(μ₂,σ₂),w₂)


# 数値積分
function approx_integration_2D(w_range, p)
    Δ = w_range[2] - w_range[1]
    (X, Y) -> sum([p(X, Y, w₁, w₂) * Δ^2 for w₁ in w_range, w₂ in w_range])
end


### シミュレーション
X_obs = [-2, 1, 5]
Y_obs = [-2.2, -1.0, 1.5]

σ = 1.0
μ₁ = 0.0
μ₂ = 0.0
σ₁ = 10.0
σ₂ = 10.0

# wの積分範囲
w_range = range(-5,5,length=500)


# 数値積分を実行
p_marginal = approx_integration_2D(w_range, p_joint)

#############################
### w₁,w₂に応じた事後分布 ###
#############################

# 事後分布の計算
w_posterior = [p_joint(X_obs, Y_obs, w₁, w₂) for w₁ in w_range, w₂ in w_range] ./ p_marginal(X_obs, Y_obs)

fig, axes = subplots(1, 2, figsize = (8, 4))

# 等高線図
cs = axes[1].contour(w_range, w_range, w_posterior', cmap = "jet")
axes[1].clabel(cs, inline = true)
set_options(axes[1], "w₁", "w₂", "posterior(contour)")

# カラーメッシュ
xgrid = repeat(w_range', length(w_range), 1)
ygrid = repeat(w_range, 1, length(w_range))
axes[2].pcolormesh(xgrid, ygrid, w_posterior', cmap = "jet", shading = "auto")
set_options(axes[2], "w₁", "w₂", "posterior(coloured)")

tight_layout()


1.1.2 予測値分布シミュレーション

 新しい入力x_pが得られたときの予測値y_pの分布を計算する。計算するのは予測分布



\begin{aligned}
P\{y_P|x_P,\boldsymbol{X},\boldsymbol{Y}\}=\displaystyle{\int P\{y_P|x_P,\boldsymbol{w}\}P\{\boldsymbol{w}|\boldsymbol{X},\boldsymbol{Y}\}}d\boldsymbol{w}
\end{aligned}


である。

######################
### 予測分布の計算 ###
######################

function approx_predictive(w_posterior, w_range, p)
    Δ = w_range[2] - w_range[1]
    (x, y) -> sum([p(x, y, w₁, w₂) * w_posterior[i, j] * Δ^2 for (i, w₁) in enumerate(w_range), (j, w₂) in enumerate(w_range)])
end

p_likelihood(xₚ,yₚ, w₁, w₂) = pdf(Normal(w₁ * xₚ + w₂, σ), yₚ)
p_predictive = approx_predictive(w_posterior, w_range, p_likelihood)

xₚ = 4.0
fig, ax = subplots()

ys = range(-5, 5, length = 100)
ax.plot(ys, p_predictive.(xₚ, ys))
set_options(ax, "yₚ", "density", "predictive distribution p(y_p|x_p=$(xₚ), X=X_obs, Y=Y_obs)")


※計算に少し時間が掛かります。

################################
### 予測値に対する密度の変化 ###
################################

# 描画範囲
xs = range(-10, 10, length = 100)
ys = range(-5, 5, length = 100)


# 密度の計算
density_y = p_predictive.(xs, ys')

fig,axes = subplots(1,2,sharey = true, figsize = (8,4))


# 等高線図
cs = axes[1].contour(xs, ys, density_y', cmap = "jet")
axes[1].clabel(cs, inline = true)
axes[1].scatter(X_obs, Y_obs)
set_options(axes[1], "x", "y", "predictive distribution(contour)")

# カラーメッシュ
xgrid = repeat(xs', length(ys), 1)
ygrid = repeat(ys, 1, length(xs))
axes[2].pcolormesh(xgrid, ygrid, w_posterior', cmap = "jet", shading = "auto")
set_options(axes[2], "w₁", "w₂", "posterior(coloured)")

tight_layout()

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