90672


現在、須山敦志「Juliaで作って学ぶベイズ統計学」(第1刷)を参照しながらJuliaでのハミルトニアン・モンテカルロ法を実装しています。ですが、P.190-195に掲載されているとおりに実装し、実例の線形回帰(同書P.196)を計算させても型異常で上手く動作しません。

残念ながらJuliaに未習熟でこれ以上自分で解決できないため、こちらで伺います。正常に動作させるためにどうしたら良いかアドバイス頂けますならば幸甚です。  
なおσ = 1.0,μ₁ = 0.0,μ₂ = 0.0,σ₁ = 10.0およびσ₂ =
10.0はエラーを起こすのが明らかだったので、本書中に記載が無いのですが自分で追加しました。また書籍の実装はVer.1.6.1で、自身はVer.1.6.2およびVer.1.8.0で検証しています。

[![画像の説明をここに入力](https://i.stack.imgur.com/LVpsj.png)](https://i.stack.imgur.com/LVpsj.png)

* * *
[code]

    using Distributions, PyPlot, ForwardDiff, LinearAlgebra
    
    # n次元単位行列
    eye(n) = Diagonal{Float64}(I, n)
    
    # グラフの諸設定を行う関数
    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
    
    # ガウス分布を提案分布としたMetropolis-Hasting法
    function GaussianMH(log_p_tilde, μ₀;
            maxiter::Int64, σ::Float64)
        # サンプルを格納する配列
        D = length(μ₀)
        μ_samples = Array{typeof(μ₀[1]), 2}(undef, D, maxiter)
    
        μ_samples[:,1] = μ₀
    
        # 受容サンプル数
        num_accepted = 1
    
        for i in 2:maxiter
        #
            μ_tmp = rand(MvNormal(μ_samples[:, i-1], σ*eye(D)))
        
            # 比率 r
            log_r = (log_p_tilde(μ_tmp) +
                     logpdf(MvNormal(μ_tmp, σ), μ_samples[:, i-1])) -
                    (log_p_tilde(μ_samples[:, i-1]) + 
                     logpdf(MvNormal(μ_samples[:,i-1], σ), μ_tmp))
            #
            is_accepted = min(1, exp(log_r)) > rand()
            new_sample = is_accepted ? μ_tmp : μ_samples[:, i-1]
        
            #
            μ_samples[:, i] = new_sample
        
            #
            num_accepted += is_accepted
        end
    
        μ_samples, num_accepted
    end
    
    # Hamiltonian Monte Carlo method
    function HMC(log_p_tilde, μ₀, maxiter::Int64, L::Int64, ε::Float64)
        # leapfrog による値の更新
        function leapflog(grad, p_in,μ_in, L, ε)
            μ=μ_in
            p = p_in + 0.5 *ε*grad(μ)
            for l in 1: L-1
                μ +=ε*p
                p += ε*grad(μ)
            end
            μ += ε*p
            p += 0.5 *ε*grad(μ)
            p, μ
        end
        
        # 非正規化対数事後分布の勾配関数
        grad(μ) = ForwardDiff.gradient(log_p_tilde, μ)
        
        # サンプルを格納する配列
        D = length(μ₀)
        μ_samples = Array{typeof(μ₀[1]), 2}(undef, D, maxiter)
        
        # 初期サンプル
        μ_samples[:, 1] = μ₀
        
        # 受容されたサンプル数
        num_accepted = 1
        
        for i in 2:maxiter
            # 運動量の生成
            p_in = randn(size(μ₀))
            
            # リープフロッグ
            p_out, μ_out = leapflog(grad, p_in, μ_samples[:, i-1], L, ε)
            
            # 比率rの対数を計算
            μ_in = μ_samples[:, i-1]
            log_r = (log_p_tilde(μ_out) + 
                     logpdf(MvNormal(zeros(D),eye(D)), vec(p_out))) -
                    (log_p_tilde(μ_in) + 
                     logpdf(MvNormal(zeros(D),eye(D)),vec(p_in)))
                    
            # 確率rでサンプル受容
            is_accepted = min(1, exp(log_r)) > rand()
            new_sample = is_accepted ? μ_out : μ_in
            
            # 新サンプルを格納
            μ_samples[:,i] = new_sample
            
            # 受容された場合、合計を加算する
            num_accepted += is_accepted
        end
        
        μ_samples, num_accepted
    end
    
    ### ラッパー関数 ###
    function inference_wrapper_GMH(log_joint, params, w_init, maxiter::Int64 = 100_000, σ::Float64 = 1.0)
        ulp(w) = log_joint(w, params)
        GaussianMH(ulp, w_init, maxiter = maxiter, σ = σ)
    end
    
    function inference_wrapper_HMC(log_joint, params, w_init,maxiter::Int64 = 100_000, L::Int64 = 100, ε::Float64 = 1e-1)
        ulp(w) = log_joint(w, params)
        HMC(ulp, w_init, maxiter = maxiter, L=L, ε=ε)
    end
    
    ###
    X_obs = [-2, 1, 5]
    Y_obs = [-2.2, -1.0, 1.5]
    
    ## 以下5つのパラメータは本書には無し
    σ = 1.0
    μ₁ = 0.0
    μ₂ = 0.0
    σ₁ = 10.0
    σ₂ = 10.0
    
    log_joint(w, X, Y, σ, μ₁,σ₁,μ₂,σ₂) =
              sum(logpdf.(Normal.(w[1]*X .+ w[2],σ), Y)) +
              logpdf(Normal(μ₁,σ₁), w[1]) +
              logpdf(Normal(μ₂,σ₂), w[2])
    params = (X_obs, Y_obs, σ, μ₁, σ₁,μ₂,σ₂)
    
    ulp(w) = log_joint(w, params...)
    
    
    
    ###
    # 初期値
    w_init = randn(2)
    
    param_posterior_GMH, num_accepted_GMH = 
        inference_wrapper_GMH(log_joint, params, w_init, maxiter = 300, σ = 1.0)
    param_posterior_HMC, num_accepted_HMC = 
        inference_wrapper_HMC(log_joint, params, w_init, maxiter = 300, L = 10, ε = 1e-1)
    
    # 
    fig, axes = subplots(2,1,figsize=(8,4))
    
    axes[1].plot(param_posterior_GMH[1,:])
    set_options(axes[1], "iteration", "w₁", "w₁ sequence (GMH)")
    
    axes[2].plot(param_posterior_GMH[2,:])
    set_options(axes[2], "iteration", "w₂", "w₂ sequence (GMH)")
    
    tight_layout()
    println("acceptance rate (GMH) = $(num_accepted_GMH/maxiter)")
    
[/code]


CC BY-SA 4.0
2022-08-21T06:38:55.140
90670
2022-08-21T07:34:41.200
54084
question
julia
ハミルトニアン・モンテカルロ法の実装関数が型異常を起こす


エラーメッセージを見ると`got unsupported keyword arguments`ですので型ではなくキーワード引数の問題のようです。

質問のコードを見ると、`inference_wrapper_GMH`のキーワード引数の前に`;`が抜けているようです。（`inference_wrapper_HMC`も同様）


CC BY-SA 4.0
2022-08-21T07:34:41.200
90672
2022-08-21T07:34:41.200
39819
90670
answer