1. 勾配を用いる必要性
対数事後分布の勾配を計算するのが今後のアプローチのコアである。
1.2 勾配を用いた計算の効率化
1.2.3 ラプラス近似
using Distributions, PyPlot, ForwardDiff, LinearAlgebra # n次元単位行列 eye(n) = Diagonal{Forward64}(I, n) # パラメータ抽出用の関数 unzip(a) = map(x -> getfield(a, x), fieldnames(eltype(a))) # グラフの諸設定を行う関数 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 ### Laplace近似 # (1) 勾配法によって事後分布の極大値を1つ求め、近似用の正規分布の平均とする # (2) 求めた極大値において2回微分が一致するように近似用の正規分布の分散を与える X_obs = [-2,1,5] Y_obs = [-2.2,-1.0,1.5] # まずは簡単にy=w₁xを推定する σ = 1.0 # 誤差項の標準偏差 # 事前分布の平均値と標準偏差 μ₁ = 0.0 σ₁ = 10.0 w₂ = 0.0 ulp(w₁) = sum(logpdf.(Normal.(w₁*X_obs .+ w₂,σ), Y_obs)) + logpdf(Normal(μ₁, σ₁), w₁) w₁s = range(-5, 5, length = 100) fig,axes = subplots(1, 2, figsize = (8,4)) # axes[1].plot(w₁s, ulp.(w₁s)) set_options(axes[1], "w₁","log density (unnormalised)", "unnormalised log posterior distribution") axes[2].plot(w₁s, exp.(ulp.(w₁s))) set_options(axes[2], "w₁","log density (unnormalised)", "unnormalised log posterior distribution") tight_layout() function gradient_method_1dim(f, x_init, η, maxiter) f′(x) = ForwardDiff.derivative(f, x) x_seq = Array{typeof(x_init), 1}(undef, maxiter) x_seq[1] = x_init for i in 2:maxiter x_seq[i] = x_seq[i-1] + η*f′(x_seq[i-1]) end x_seq end # 最適化パラメータ w₁_init = 0.0 maxiter = 100 η = 0.01 # 最適化の実施 w₁_seq = gradient_method_1dim(ulp, w₁_init, η, maxiter) # 勾配法の過程を可視化 fig, ax = subplots(figsize = (8,4)) ax.plot(w₁_seq) set_options(ax, "iteration","w₁", "w₁ sequence") ### ラッパー関数を挟む ### :関数外の変数を直接参照するとJuliaの計算速度が著しく低下する ### :そこでラッパー関数を挟んで、計算に必要なデータやパラメータをすべて関数の引数として与える # 最適化パラメータ w₁_init = 0.0 maxiter = 100 η = 0.01 # 最適化ラッパー関数の定義 function inference_wrapper_gd_1dim(log_joint, params, w₁_init, η, maxiter) ulp(w₁) = log_joint(w₁, params...) w₁_seq = gradient_method_1dim(ulp, w₁_init, η, maxiter) w₁_seq end # 対数同時分布 log_joint(w₁, X, Y, w₂,σ,μ₁,σ₁) = sum(logpdf.(Normal.(w₁*X.+w₂,σ), Y)) + logpdf(Normal(μ₁,σ₁), w₁) params = (X_obs, Y_obs, w₁_init, σ, μ₁, σ₁) # 最適化の実施 w₁_seq = inference_wrapper_gd_1dim(log_joint, params, w₁_init, η, maxiter) ############################################################################### ### 求めた極大値において2階微分が一致するように近似用正規分布の分散を求める ### ############################################################################### # 近似分布用の平均 μ_approx = w₁_seq[end] # 近似用正規分布の分散 grad(x) = ForwardDiff.derivative(ulp, x) hessian(x) = ForwardDiff.derivative(grad, x) σ_approx = sqrt(inv(-hessian(μ_approx))) # 近似用正規分布を図示してみる w₁s = range(-0.5,1,length = 1000) fig,axes = subplots(1,2,figsize = (8,4)) axes[1].plot(w₁s, exp.(ulp.(w₁s))) set_options(axes[1], "w₁", "density (unnormalised)", "unnormalised posterioir distribution") # 得られた近似分布の可視化 axes[2].plot(w₁s, pdf.(Normal.(μ_approx, σ_approx), w₁s)) set_options(axes[2], "w₁", "density", "approximate distribution") tight_layout()