前回
1. 数値積分を使った線形モデルの学習
1.1 理論の準備
条件付き分布の定義から、学習用のデータセットおよび
が与えられた後の母数の事後分布は
で与えることができる。ここで
である。したがって上記の積分計算ができるならば、事後分布の密度関数が計算できる。
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 予測値分布シミュレーション
新しい入力が得られたときの予測値
の分布を計算する。計算するのは予測分布
である。
###################### ### 予測分布の計算 ### ###################### 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()
