Juliaを使ってみる(16/22)：ラプラス近似(その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()

次回

power-of-awareness.com

今回やりたいこと

前回

1. 勾配を用いる必要性

1.2 勾配を用いた計算の効率化

1.2.3 ラプラス近似

次回

1.　勾配を用いる必要性

1.2　勾配を用いた計算の効率化

1.2.3　ラプラス近似