りんごがでている

何か役に立つことを書きます

JuliaTokyo #2 でモンテカルロ法について話してきました

Julia

9月27日の JuliaTokyo #2 で『Juliaで学ぶ Hamiltonian Monte Carlo (NUTS 入り)』と題してHamiltonian Monte Carlo (HMC) とStanで使われているNUTSについてお話してきました。

JuliaTokyo #2 - connpass

内容

サンプルコード: bicycle1885/JuliaTokyo2HMC · GitHub

MCMC自体については30分ではほとんど説明できませんでしたが、会場のほとんどの方が既に知っているようでしたのでちょうど良かったようです。

全体の流れとしては、Metropolis-Hastings、Hamiltonian Monte Carlo、No-U-Turn Samplerの順に3つサンプラーを実装コードと正規分布からのサンプリングを例に紹介し、どのように前のサンプラーの弱点を克服していってるかを説明しました。

Unicodeプログラミング

サンプルコードはJuliaの愉しみのひとつであるUnicode文字をふんだんに使っていたところが面白いところです。

function build_tree(L::Function, ∇L::Function, θ::Vector{Float64}, r::Vector{Float64}, u::Float64, v::Int, j::Int, ϵ::Float64)
    if j == 0
        θ′, r′ = leapfrog(∇L, θ, r, v * ϵ)
        C′ = u ≤ exp(L(θ′) - r′ ⋅ r′ / 2) ? Set([(θ′, r′)]) : Set([])
        s′ = int(L(θ′) - r′ ⋅ r′ / 2 > log(u) - Δmax)
        return θ′, r′, θ′, r′, C′, s′
    else
        θ⁻, r⁻, θ⁺, r⁺, C′, s′ = build_tree(L, ∇L, θ, r, u, v, j - 1, ϵ)
        if v == -1
            θ⁻, r⁻, _, _, C″, s″ = build_tree(L, ∇L, θ⁻, r⁻, u, v, j - 1, ϵ)
        else
            _, _, θ⁺, r⁺, C″, s″ = build_tree(L, ∇L, θ⁺, r⁺, u, v, j - 1, ϵ)
        end
        s′ = s′ * s″ * ((θ⁺ - θ⁻) ⋅ r⁻ ≥ 0) * ((θ⁺ - θ⁻) ⋅ r⁺ ≥ 0)
        C′ = C′ ∪ C″
        return θ⁻, r⁻, θ⁺, r⁺, C′, s′
    end
end

箴言

+α

社内で改造しているようですが、手続き上のアレでプルリクは受けられないようです(´・ω・`)