Turing Guide Basic Usage¶

In [1]:
using Turing
using StatsPlots
using Random
using Suppressor
using DataFrames
In [2]:
Random.seed!(1234)
Out[2]:
TaskLocalRNG()
In [3]:
@model function gdemo(x::AbstractArray{<:Real})
    # Prior of Target Parameter
    # P(s²,m) = P(s²)P(m)
    # P(s²) ∝ pdf(InverseGamma(2,3),s²)
    s² ~ InverseGamma(2,3)
    # P(m) ∝ pdf(Normal(0, √s²),s²)
    m ~ Normal(0, √s²)
    
    # Likelihood
    # x = {x₁,x₂,...,xₙ}    
    # P(x|m,s²)
    # P(x|m,s²) ∝ pdf.(Normal(m, √s²),x)    
    x ~ Normal(m, √s²)
        
    #return
    # pdf.(Normal(m, √s²),x)*Normal(0, √s²)*pdf(InverseGamma(2,3),s²)
#     return x
end
Out[3]:
gdemo (generic function with 2 methods)
In [4]:
# Make test sample
# x ~ Normal(3.5, 2.0)
x = rand(Normal(2.5, √2.0), 300)
density(x, label="Probability Density of x")
Out[4]:
In [5]:
samples = sample(gdemo(x), NUTS(),1000)

┌ Info: Found initial step size
└   ϵ = 0.2
Sampling: 100%|█████████████████████████████████████████| Time: 0:00:00
Out[5]:
Chains MCMC chain (1000×14×1 Array{Float64, 3}):

Iterations        = 501:1:1500
Number of chains  = 1
Samples per chain = 1000
Wall duration     = 10.03 seconds
Compute duration  = 10.03 seconds
parameters        = s², m
internals         = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size

Summary Statistics
  parameters      mean       std      mcse    ess_bulk   ess_tail      rhat    ⋯
      Symbol   Float64   Float64   Float64     Float64    Float64   Float64    ⋯

          s²    2.3377    0.1950    0.0065    955.6722   640.0229    1.0081    ⋯
           m    2.3511    0.0888    0.0027   1061.0144   843.2654    0.9995    ⋯
                                                                1 column omitted

Quantiles
  parameters      2.5%     25.0%     50.0%     75.0%     97.5% 
      Symbol   Float64   Float64   Float64   Float64   Float64 

          s²    2.0204    2.2034    2.3190    2.4602    2.7617
           m    2.1845    2.2896    2.3513    2.4093    2.5246
In [6]:
fieldnames(typeof(samples))
Out[6]:
(:value, :logevidence, :name_map, :info)
In [7]:
samples.name_map
Out[7]:
(parameters = [:s², :m], internals = [:lp, :n_steps, :is_accept, :acceptance_rate, :log_density, :hamiltonian_energy, :hamiltonian_energy_error, :max_hamiltonian_energy_error, :tree_depth, :numerical_error, :step_size, :nom_step_size])
In [8]:
plot(samples)
Out[8]:
In [9]:
samples[:m][1:3] |> display

1-dimensional AxisArray{Float64,1,...} with axes:
    :row, Base.OneTo(3)
And data, a 3-element Vector{Float64}:
 2.5001000822783443
 2.3754059003239014
 2.355634943865946
In [10]:
samples[:s²][1:3] |> display
samples["s²"][1:3]|> display

1-dimensional AxisArray{Float64,1,...} with axes:
    :row, Base.OneTo(3)
And data, a 3-element Vector{Float64}:
 2.374601897723656
 2.1637177928075286
 2.2423605153574666
1-dimensional AxisArray{Float64,1,...} with axes:
    :row, Base.OneTo(3)
And data, a 3-element Vector{Float64}:
 2.374601897723656
 2.1637177928075286
 2.2423605153574666

Sampling Multiple Chains¶

Multithreaded sampling¶

In [11]:
# sample(model, sampler, parallel_type, n, n_chains)
samples = sample(gdemo(x),NUTS(),MCMCThreads(),1000,5)
p1 = plot(samples)
fieldnames(typeof(samples)) |> display
size(samples) |> display
size(samples[:m]) |> display
(mean(samples[:m]) , mean(samples[:s²])) |> display
plot(p1)

┌ Info: Found initial step size
└   ϵ = 0.05
┌ Info: Found initial step size
└   ϵ = 0.4
┌ Info: Found initial step size
└   ϵ = 0.2
┌ Info: Found initial step size
└   ϵ = 0.05
┌ Info: Found initial step size
└   ϵ = 0.0125

(:value, :logevidence, :name_map, :info)
(1000, 14, 5)
(1000, 5)
(2.3540283799147232, 2.3274197580434497)
Out[11]:
In [12]:
samples[:,:,1]
Out[12]:
Chains MCMC chain (1000×14×1 Array{Float64, 3}):

Iterations        = 501:1:1500
Number of chains  = 1
Samples per chain = 1000
Wall duration     = 0.48 seconds
Compute duration  = 1.28 seconds
parameters        = s², m
internals         = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size

Summary Statistics
  parameters      mean       std      mcse    ess_bulk   ess_tail      rhat    ⋯
      Symbol   Float64   Float64   Float64     Float64    Float64   Float64    ⋯

          s²    2.3299    0.1875    0.0062    932.9716   654.7551    1.0000    ⋯
           m    2.3572    0.0896    0.0026   1187.4449   735.2058    1.0001    ⋯
                                                                1 column omitted

Quantiles
  parameters      2.5%     25.0%     50.0%     75.0%     97.5% 
      Symbol   Float64   Float64   Float64   Float64   Float64 

          s²    1.9799    2.2008    2.3186    2.4470    2.7127
           m    2.1927    2.2958    2.3554    2.4200    2.5372
In [13]:
@model function gdemo(x)
    n = length(x)
    s²= Vector{Real}(undef,n)
    m  = Vector{Real}(undef,n)
    for i in eachindex(x)            
        s²[i] ~ InverseGamma(2,3)
        m[i] ~ Normal(0, √(s²[i]))            
        x[i] ~ Normal(m[i],√(s²[i]))
    end
end

x = rand(Normal(2.5, √2.0), 200)
y = rand(Normal(-1.0, √3.0), 200)
p1 = density(x, label="Probability Density of x")
density!(p1,y, label="Probability Density of y") |> display


model = gdemo([x,y])

samples = sample(model, NUTS(), MCMCThreads(), 1000, 5)
samples |> display
plot(samples)

┌ Info: Found initial step size
└   ϵ = 0.05
┌ Info: Found initial step size
└   ϵ = 0.05
┌ Info: Found initial step size
└   ϵ = 0.05
┌ Info: Found initial step size
└   ϵ = 0.2
┌ Info: Found initial step size
└   ϵ = 0.05

Chains MCMC chain (1000×16×5 Array{Float64, 3}):

Iterations        = 501:1:1500
Number of chains  = 5
Samples per chain = 1000
Wall duration     = 5.01 seconds
Compute duration  = 24.59 seconds
parameters        = s²[1], s²[2], m[1], m[2]
internals         = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size

Summary Statistics
  parameters      mean       std      mcse    ess_bulk    ess_tail      rhat   ⋯
      Symbol   Float64   Float64   Float64     Float64     Float64   Float64   ⋯

       s²[1]    1.6677    0.1632    0.0020   6988.8789   4163.0998    1.0007   ⋯
       s²[2]    3.3253    0.3276    0.0039   7113.3700   4380.5059    1.0003   ⋯
        m[1]    2.2191    0.0920    0.0010   7790.1820   4024.6893    1.0005   ⋯
        m[2]   -0.8220    0.1277    0.0015   7296.1149   3641.5370    0.9999   ⋯
                                                                1 column omitted

Quantiles
  parameters      2.5%     25.0%     50.0%     75.0%     97.5% 
      Symbol   Float64   Float64   Float64   Float64   Float64 

       s²[1]    1.3851    1.5530    1.6569    1.7733    2.0131
       s²[2]    2.7434    3.0999    3.3094    3.5232    4.0239
        m[1]    2.0405    2.1579    2.2176    2.2807    2.4021
        m[2]   -1.0749   -0.9065   -0.8207   -0.7358   -0.5739
Out[13]:
In [14]:
samples["m[1]"][1:2] |> display
group(samples,"m") |> display

1-dimensional AxisArray{Float64,1,...} with axes:
    :row, Base.OneTo(2)
And data, a 2-element Vector{Float64}:
 2.150920623843303
 2.2805898485444533
Chains MCMC chain (1000×2×5 Array{Float64, 3}):

Iterations        = 501:1:1500
Number of chains  = 5
Samples per chain = 1000
Wall duration     = 5.01 seconds
Compute duration  = 24.59 seconds
parameters        = m[1], m[2]
internals         = 

Summary Statistics
  parameters      mean       std      mcse    ess_bulk    ess_tail      rhat   ⋯
      Symbol   Float64   Float64   Float64     Float64     Float64   Float64   ⋯

        m[1]    2.2191    0.0920    0.0010   7790.1820   4024.6893    1.0005   ⋯
        m[2]   -0.8220    0.1277    0.0015   7296.1149   3641.5370    0.9999   ⋯
                                                                1 column omitted

Quantiles
  parameters      2.5%     25.0%     50.0%     75.0%     97.5% 
      Symbol   Float64   Float64   Float64   Float64   Float64 

        m[1]    2.0405    2.1579    2.2176    2.2807    2.4021
        m[2]   -1.0749   -0.9065   -0.8207   -0.7358   -0.5739

Distributed sampling (Multiple processes)¶

In [15]:
using Distributed

let    
    n_process = 4
    addprocs(n_process)
    
    # 모든 process에서 turing과 model을 사용할 수 있게
    @everywhere  using Turing
    
    @everywhere  @model function gdemo(x)
        n = length(x)
        s²= Vector{Real}(undef,n)
        m  = Vector{Real}(undef,n)
        for i in eachindex(x)            
            s²[i] ~ InverseGamma(2,3)
            m[i] ~ Normal(0, √(s²[i]))            
            x[i] ~ Normal(m[i],√(s²[i]))
        end
    end

    @everywhere  x = rand(Normal(2.5, √2.0), 200)
    @everywhere  y = rand(Normal(-1.0, √3.0), 200)
    p1 = density(x, label="Probability Density of x")
    density!(p1,y, label="Probability Density of y") |> display


    @everywhere  model = gdemo([x,y])

    samples = sample(model, NUTS(), MCMCThreads(), 1000, 5)
    samples |> display
    plot(samples)    
end

┌ Info: Found initial step size
└   ϵ = 0.025
┌ Info: Found initial step size
└   ϵ = 0.025
┌ Info: Found initial step size
└   ϵ = 0.05
┌ Info: Found initial step size
└   ϵ = 0.025
┌ Info: Found initial step size
└   ϵ = 0.05

Chains MCMC chain (1000×16×5 Array{Float64, 3}):

Iterations        = 501:1:1500
Number of chains  = 5
Samples per chain = 1000
Wall duration     = 1.4 seconds
Compute duration  = 5.67 seconds
parameters        = s²[1], s²[2], m[1], m[2]
internals         = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size

Summary Statistics
  parameters      mean       std      mcse    ess_bulk    ess_tail      rhat   ⋯
      Symbol   Float64   Float64   Float64     Float64     Float64   Float64   ⋯

       s²[1]    1.8682    0.1821    0.0021   7532.7276   3863.6923    1.0005   ⋯
       s²[2]    2.6048    0.2585    0.0028   8597.8616   4103.2110    0.9999   ⋯
        m[1]    2.5091    0.0954    0.0010   9175.4662   3751.2857    1.0003   ⋯
        m[2]   -0.9552    0.1111    0.0013   7586.4452   4357.6853    1.0018   ⋯
                                                                1 column omitted

Quantiles
  parameters      2.5%     25.0%     50.0%     75.0%     97.5% 
      Symbol   Float64   Float64   Float64   Float64   Float64 

       s²[1]    1.5476    1.7375    1.8564    1.9860    2.2552
       s²[2]    2.1437    2.4223    2.5874    2.7659    3.1640
        m[1]    2.3279    2.4458    2.5090    2.5723    2.6942
        m[2]   -1.1709   -1.0285   -0.9553   -0.8812   -0.7372
Out[15]:

Sampling from an Unconditional Distribution(The Prior) - 사전분포에서 샘플링 하기¶

  • 사전분포로 부터 샘플링 하기
  • 관측값(x) 을 missing으로 넘기면 사전분포로 부터 샘플링
  • 입력이나 샘플러를 지정하지 않고 모델을 호출하여 사전분포에서 모델을 (마치 함수처럼) 실행할 수도 있습니다. 아래 예제에서는 x와 y라는 두 변수를 반환하는 gdemo 모델을 지정합니다. 이 모델에는 x와 y가 인수로 포함되지만, x나 y를 전달하지 않고 함수를 호출하면 Turing의 컴파일러는 해당 분포(사후분포)에서 가져올 값이 누락된 것으로 간주합니다. 반환 문은 사전분포로 부터 샘플링된 x 및 y 값을 검색하는 데 필요합니다.
  • ※주의 : 모델에서 사용한 파라미터 m, s²를 vector로 사용하여 push!(m,v[1]) 을 호출하면 에러가 발생 하는데 이유는 모르겠으나 문제 발생하니 model에서 사용한 파라미터(모수)는 파해야 함

$ P(m,s^2,x_1,x_2) \propto P(x_1) \cdot P(x_2) \cdot P(m) \cdot P(s^2)\\ $

In [16]:
using LaTeXStrings
let
    @model function gdemo(x,y)
        s² ~ InverseGamma(2,3)
        m ~ Normal(0, √s²)
        x ~ Normal(m, √s²)
        y ~ Normal(m, √s²)
        return x,y
    end
    
    # model을 함수 처럼 직접 호출하여 사전분포에서 sampling
    # ※주의 : 모델에서 사용한 파라미터 m, s²를 vector로 
    #          사용하여 push!(m,v[1]) 을 호출하면 에러가 발생
    #          하는데 이유는 모르겠으나 문제 발생하니
    #          model에서 사용한 파라미터(모수)는 파해야 함
    g_prior_sample = gdemo(missing, missing)
    σ² = Vector{Real}()
    μ= Vector{Real}()
    
    n_sample  = 1000
    
    for i in 1:n_sample
        v = g_prior_sample()
        push!(σ²,v[1])
        push!(μ,v[2])
    end    
    p1 = density(σ²,label=L"s^2")
    p2 = density(μ, label="m")
    "model을 함수처럼 호출하여사전분포에서 샘플링 하기" |> display
    plot(p1,p2,layout=(1,2)) |> display
    
    # sampler 통해 사전분포에서 샘플링 하기
    model = gdemo(missing, missing)
    samples = sample(model, Prior(), n_sample)
    "Sampler를 통해 사전분포에서 샘플링 하기" |> display
    plot(samples) |> display
    samples
end

"model을 함수처럼 호출하여사전분포에서 샘플링 하기"
"Sampler를 통해 사전분포에서 샘플링 하기"
Out[16]:
Chains MCMC chain (1000×5×1 Array{Float64, 3}):

Iterations        = 1:1:1000
Number of chains  = 1
Samples per chain = 1000
Wall duration     = 0.31 seconds
Compute duration  = 0.31 seconds
parameters        = s², m, x, y
internals         = lp

Summary Statistics
  parameters      mean       std      mcse    ess_bulk    ess_tail      rhat   ⋯
      Symbol   Float64   Float64   Float64     Float64     Float64   Float64   ⋯

          s²    3.1413    6.7636    0.2198    893.6556    907.7013    1.0011   ⋯
           m    0.0544    1.6973    0.0528   1008.9203    889.9348    0.9999   ⋯
           x   -0.0334    2.5532    0.0777   1080.0696   1030.6930    0.9995   ⋯
           y    0.0303    2.3488    0.0737    985.3144   1025.7469    0.9994   ⋯
                                                                1 column omitted

Quantiles
  parameters      2.5%     25.0%     50.0%     75.0%     97.5% 
      Symbol   Float64   Float64   Float64   Float64   Float64 

          s²    0.5659    1.1395    1.7994    3.3211   13.1305
           m   -3.5274   -0.8647    0.0864    0.9247    3.7447
           x   -5.3766   -1.4937    0.0084    1.2415    5.3605
           y   -4.5638   -1.2242   -0.0006    1.3849    4.9284

Sampling from a Conditional Distribution (The Posterior) - 파라미터를 사후분포에서 샘플링 하기¶

Treating observations as random variables¶

  • 모델에 missing 값이 있는 입력값은 추정/샘플링할 파라미터(모수), 즉 랜덤 변수(Random Variables)로 취급됩니다. 이는 해당 파라미터에 대한 추출을 시뮬레이션하거나 조건부 분포에서 샘플링하려는 경우에 유용할 수 있습니다. 튜링은 다음 구문을 지원합니다.
  • 이런 경우 파라미터를 사후분포에서 샘플링함
In [17]:
let
    @model function gdemo(x, ::Type{T} = Float64) where {T}
        # 이 경우 x는 데이터가 아니라 파라미터로 취급 된다.
        if x === missing
            # Initialize `x` if missing
            x = Vector{T}(undef,2)
        end
        s²~ InverseGamma(2,3)
        m ~ Normal(0, √s²)
        for i in eachindex(x)
            x[i] ~ Normal(m, √s²)
        end
    end
    
    # Construct a model with x = missing
    model = gdemo(missing)
#     c = sample(model, HMC(0.01,5), 500)    
    c = sample(model, NUTS(), 500)    
    plot(c) |> display
    c |> display    
end

┌ Info: Found initial step size
└   ϵ = 0.42500000000000004
Sampling: 100%|█████████████████████████████████████████| Time: 0:00:00

Chains MCMC chain (500×16×1 Array{Float64, 3}):

Iterations        = 251:1:750
Number of chains  = 1
Samples per chain = 500
Wall duration     = 4.15 seconds
Compute duration  = 4.15 seconds
parameters        = s², m, x[1], x[2]
internals         = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size

Summary Statistics
  parameters      mean       std      mcse   ess_bulk   ess_tail      rhat   e ⋯
      Symbol   Float64   Float64   Float64    Float64    Float64   Float64     ⋯

          s²    2.7624    2.9191    0.2679    99.8749    94.0338    0.9982     ⋯
           m   -0.0045    1.5138    0.1159   178.1741   172.4582    1.0088     ⋯
        x[1]   -0.0850    2.2217    0.1396   267.7089   210.8655    1.0056     ⋯
        x[2]    0.0700    2.1855    0.1436   231.3029   168.2802    1.0030     ⋯
                                                                1 column omitted

Quantiles
  parameters      2.5%     25.0%     50.0%     75.0%     97.5% 
      Symbol   Float64   Float64   Float64   Float64   Float64 

          s²    0.5052    1.1032    1.9166    3.1812   10.9506
           m   -3.7048   -0.8142   -0.0125    0.8082    3.4494
        x[1]   -5.0158   -1.2540   -0.0730    1.1185    4.3958
        x[2]   -4.5940   -1.1276    0.1646    1.1763    4.4872
또한 튜링은 x에서 `missing` 값과 non-`missing` 값이 혼합된 경우, `missing`은 랜덤 변수로 처리하여 샘플링하고 나머지 값은 관측값으로 처리하는 기능도 지원합니다.
$ P(m,s^2,x_1|x_2) \propto P(x_2|m,s^2,x_1) \cdot P(m,s^2,x_1) \\ \propto P(x_2|m) \cdot P(x_2|s^2) \cdot P(x_2|x_1) \cdot P(m) \cdot P(s^2) \cdot P(x_1) \\ $
In [18]:
let
    @model function gdemo(x)
        s²~ InverseGamma(2, 3)
        m ~ Normal(0, √s²)
        for i in eachindex(x)
            x[i] ~ Normal(m, √s²)
        end
    end
    # x[1] is parameter, but x[2] is an observation
    model = gdemo([missing, 2.4])
#     c = sample(model, HMC(0.01, 5), 1000) 
    c = sample(model, NUTS(), 1000) 
    plot(c) |> display
    c |> display
end

┌ Info: Found initial step size
└   ϵ = 0.8
Sampling: 100%|█████████████████████████████████████████| Time: 0:00:00

Chains MCMC chain (1000×15×1 Array{Float64, 3}):

Iterations        = 501:1:1500
Number of chains  = 1
Samples per chain = 1000
Wall duration     = 3.8 seconds
Compute duration  = 3.8 seconds
parameters        = s², m, x[1]
internals         = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size

Summary Statistics
  parameters      mean       std      mcse   ess_bulk   ess_tail      rhat   e ⋯
      Symbol   Float64   Float64   Float64    Float64    Float64   Float64     ⋯

          s²    2.6711    2.3203    0.1163   343.5518   378.3130    0.9998     ⋯
           m    1.2877    1.1422    0.0564   429.2675   412.8486    1.0017     ⋯
        x[1]    1.4401    1.9348    0.1014   376.5340   393.0771    1.0009     ⋯
                                                                1 column omitted

Quantiles
  parameters      2.5%     25.0%     50.0%     75.0%     97.5% 
      Symbol   Float64   Float64   Float64   Float64   Float64 

          s²    0.7772    1.2832    1.9731    3.1616    9.1785
           m   -0.9917    0.6002    1.3167    1.9541    3.6421
        x[1]   -2.3558    0.2496    1.3502    2.5362    5.5510

Default Values¶

튜링 모델에 대한 인수는 일반 줄리아 함수에서 기본값이 작동하는 방식과 매우 유사하게 기본값을 가질 수 있습니다. 예를 들어 다음은 x에 missing을 할당하고 이를 랜덤 변수로 처리합니다. 기본값이 누락되지 않은 경우 x에 해당 값이 할당되고 대신 관측값으로 처리됩니다.

In [19]:
let
    @model function generative(x = missing, ::Type{T} = Float64) where {T <: Real}
        if x === missing
            # Initialize x when missing
            x = Vector{T}(undef,10)
        end
        s² ~ InverseGamma(2, 3)
        m ~ Normal(0, √s²)
        for i in 1:length(x)
            x[i] ~ Normal(m, √s²)
        end
        return s²,m
    end
    
    model = generative()
    samples = sample(model, NUTS(), 1000)
    plot(samples)
end

┌ Info: Found initial step size
└   ϵ = 0.8
Sampling: 100%|█████████████████████████████████████████| Time: 0:00:00
Out[19]:

Access Values inside Chain¶

체인 내부의 값에 액세스하는 방법은 여러 가지가 있습니다:

  1. DataFrame 개체로 변환
  2. raw AxisArray 사용
  3. 3차원 Array 객체 만들기

예를 들어, c가 체인이라고 가정합니다:

  1. DataFrame(c)는 c를 DataFrame으로 변환합니다,
  2. c.value는 c 내부의 값을 AxisArray로, 그리고
  3. c.value.data는 c 내부의 값을 3D Array로 검색합니다.
In [20]:
@model function gdemo(x)
    n = length(x)
    s²= Vector{Real}(undef,n)
    m  = Vector{Real}(undef,n)
    for i in eachindex(x)            
        s²[i] ~ InverseGamma(2,3)
        m[i] ~ Normal(0, √(s²[i]))            
        x[i] ~ Normal(m[i],√(s²[i]))
    end
end

x = rand(Normal(2.5, √2.0), 200)
y = rand(Normal(-1.0, √3.0), 200)

model = gdemo([x,y])

c = sample(model, NUTS(), MCMCThreads(), 1000, 2)
c |> display
plot(c) |> display
autocorplot(c) |> display

┌ Info: Found initial step size
└   ϵ = 0.2
┌ Info: Found initial step size
└   ϵ = 0.2

Chains MCMC chain (1000×16×2 Array{Float64, 3}):

Iterations        = 501:1:1500
Number of chains  = 2
Samples per chain = 1000
Wall duration     = 1.21 seconds
Compute duration  = 2.36 seconds
parameters        = s²[1], s²[2], m[1], m[2]
internals         = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size

Summary Statistics
  parameters      mean       std      mcse    ess_bulk    ess_tail      rhat   ⋯
      Symbol   Float64   Float64   Float64     Float64     Float64   Float64   ⋯

       s²[1]    2.0265    0.2043    0.0039   2848.5842   1573.2193    1.0000   ⋯
       s²[2]    3.0390    0.2976    0.0055   2901.4464   1667.9545    1.0004   ⋯
        m[1]    2.4018    0.1014    0.0020   2555.9499   1485.1916    1.0014   ⋯
        m[2]   -0.9410    0.1228    0.0020   3698.3345   1534.9285    0.9998   ⋯
                                                                1 column omitted

Quantiles
  parameters      2.5%     25.0%     50.0%     75.0%     97.5% 
      Symbol   Float64   Float64   Float64   Float64   Float64 

       s²[1]    1.6669    1.8835    2.0114    2.1517    2.4632
       s²[2]    2.5013    2.8319    3.0268    3.2217    3.6829
        m[1]    2.2116    2.3310    2.4014    2.4714    2.5981
        m[2]   -1.1769   -1.0262   -0.9389   -0.8613   -0.6992
In [21]:
let
    using DataFrames
    
    df = DataFrame(c)
    "DataFrame" |> display
    df |> display
    
    "c.value" |> display
    c.value |> display
    
    "c.value.data" |> display
    c.value.data |> display
end

"DataFrame"
2000×18 DataFrame
1975 rows omitted
Rowiterationchains²[1]s²[2]m[1]m[2]lpn_stepsis_acceptacceptance_ratelog_densityhamiltonian_energyhamiltonian_energy_errormax_hamiltonian_energy_errortree_depthnumerical_errorstep_sizenom_step_size
Int64Int64Float64Float64Float64Float64Float64Float64Float64Float64Float64Float64Float64Float64Float64Float64Float64Float64
150111.791712.92592.47615-1.00088-754.7833.01.00.989894-754.783755.10.0307885-0.1708082.00.00.9437740.943774
250211.791712.92592.47615-1.00088-754.7833.01.00.507124-754.783758.4730.00.7748262.00.00.9437740.943774
350311.783653.201172.44422-0.715697-756.2683.01.00.587669-756.268758.10.4347440.6386662.00.00.9437740.943774
450412.041882.718772.36068-0.952244-754.3227.01.00.971073-754.322756.786-0.53635-0.536352.00.00.9437740.943774
550512.00083.348442.48828-0.923913-754.6483.01.00.873075-754.648755.730.06956490.2746972.00.00.9437740.943774
650612.044232.641082.39721-0.955044-754.5957.01.00.624175-754.595758.041-0.03852710.9540112.00.00.9437740.943774
750712.020323.600372.29273-0.911933-755.933.01.00.347244-755.93761.5170.5499191.534662.00.00.9437740.943774
850812.110522.624962.50005-0.97364-755.2793.01.00.861626-755.279757.751-0.004677440.536352.00.00.9437740.943774
950911.963773.466452.45473-1.09452-755.5457.01.00.867588-755.545757.1310.1017350.2698722.00.00.9437740.943774
1051012.298623.6442.55443-0.889306-757.5153.01.00.746443-757.515759.2550.4077490.5548072.00.00.9437740.943774
1151111.769883.49022.38863-0.820608-756.0213.01.01.0-756.021758.849-0.282026-0.4045012.00.00.9437740.943774
1251211.948393.016242.55037-0.704091-756.697.01.00.768924-756.69759.3320.2258510.7491282.00.00.9437740.943774
1351312.077742.849752.38127-0.993518-754.0253.01.00.87768-754.025758.076-0.67024-0.670241.00.00.9437740.943774
⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮
1989148922.053963.580512.13081-1.01232-759.0233.01.01.0-759.023760.536-0.00138792-0.7378952.00.00.9635390.963539
1990149022.102313.032912.20647-1.07199-756.2553.01.01.0-756.255759.247-0.682672-0.8776642.00.00.9635390.963539
1991149122.102313.032912.20647-1.07199-756.2551.01.00.652801-756.255757.8430.00.4264831.00.00.9635390.963539
1992149221.997413.060722.59206-0.80927-756.0553.01.01.0-756.055757.351-0.0692518-0.2916392.00.00.9635390.963539
1993149321.546513.04222.37786-1.03687-757.7443.01.00.557181-757.744761.694-0.04707921.450612.00.00.9635390.963539
1994149422.164043.12452.53922-0.777052-755.793.01.00.984179-755.79759.128-0.315962-0.3159622.00.00.9635390.963539
1995149522.061722.9462.29504-1.04885-754.7313.01.00.755407-754.731759.38-0.4283270.7288622.00.00.9635390.963539
1996149621.920623.12042.52605-0.775365-755.5073.01.00.711448-755.507757.5490.3231410.8904492.00.00.9635390.963539
1997149722.141083.028752.51964-1.02001-754.7683.01.00.889048-754.768756.74-0.2394460.239822.00.00.9635390.963539
1998149821.946672.990992.28703-0.833088-754.8143.01.00.91457-754.814755.9640.05575820.225722.00.00.9635390.963539
1999149922.085832.85712.53715-0.91628-754.7813.01.00.989534-754.781755.584-0.01363380.03190162.00.00.9635390.963539
2000150022.086593.724132.3523-0.85725-756.3193.01.00.417331-756.319759.250.4362381.24052.00.00.9635390.963539
"c.value"
3-dimensional AxisArray{Float64,3,...} with axes:
    :iter, 501:1:1500
    :var, [Symbol("s²[1]"), Symbol("s²[2]"), Symbol("m[1]"), Symbol("m[2]"), :lp, :n_steps, :is_accept, :acceptance_rate, :log_density, :hamiltonian_energy, :hamiltonian_energy_error, :max_hamiltonian_energy_error, :tree_depth, :numerical_error, :step_size, :nom_step_size]
    :chain, 1:2
And data, a 1000×16×2 Array{Float64, 3}:
[:, :, 1] =
 1.79171  2.9259   2.47615  -1.00088   …  2.0  0.0  0.943774  0.943774
 1.79171  2.9259   2.47615  -1.00088      2.0  0.0  0.943774  0.943774
 1.78365  3.20117  2.44422  -0.715697     2.0  0.0  0.943774  0.943774
 2.04188  2.71877  2.36068  -0.952244     2.0  0.0  0.943774  0.943774
 2.0008   3.34844  2.48828  -0.923913     2.0  0.0  0.943774  0.943774
 2.04423  2.64108  2.39721  -0.955044  …  2.0  0.0  0.943774  0.943774
 2.02032  3.60037  2.29273  -0.911933     2.0  0.0  0.943774  0.943774
 2.11052  2.62496  2.50005  -0.97364      2.0  0.0  0.943774  0.943774
 1.96377  3.46645  2.45473  -1.09452      2.0  0.0  0.943774  0.943774
 2.29862  3.644    2.55443  -0.889306     2.0  0.0  0.943774  0.943774
 1.76988  3.4902   2.38863  -0.820608  …  2.0  0.0  0.943774  0.943774
 1.94839  3.01624  2.55037  -0.704091     2.0  0.0  0.943774  0.943774
 2.07774  2.84975  2.38127  -0.993518     1.0  0.0  0.943774  0.943774
 ⋮                                     ⋱                      ⋮
 1.84218  3.05626  2.29374  -0.911148     2.0  0.0  0.943774  0.943774
 1.99603  3.36006  2.40845  -1.07951      2.0  0.0  0.943774  0.943774
 1.98988  2.53361  2.374    -0.787454  …  3.0  0.0  0.943774  0.943774
 1.98688  2.67816  2.35937  -1.17963      2.0  0.0  0.943774  0.943774
 1.74054  2.76433  2.4468   -0.91494      2.0  0.0  0.943774  0.943774
 1.6844   3.23351  2.39102  -0.881826     2.0  0.0  0.943774  0.943774
 1.70373  2.93048  2.36806  -1.00723      2.0  0.0  0.943774  0.943774
 1.95558  3.21013  2.29591  -0.785725  …  2.0  0.0  0.943774  0.943774
 2.25762  3.5584   2.14025  -0.684222     2.0  0.0  0.943774  0.943774
 2.02272  3.16815  2.28926  -0.661574     2.0  0.0  0.943774  0.943774
 1.98888  2.85139  2.28506  -1.10707      2.0  0.0  0.943774  0.943774
 2.22531  2.53992  2.42093  -0.805889     2.0  0.0  0.943774  0.943774

[:, :, 2] =
 1.72913  3.374    2.48625  -0.850005  …  2.0  0.0  0.963539  0.963539
 2.34423  2.83456  2.1521   -0.915034     2.0  0.0  0.963539  0.963539
 2.12179  3.28248  2.55795  -0.98139      2.0  0.0  0.963539  0.963539
 1.93434  2.76009  2.25582  -0.898026     2.0  0.0  0.963539  0.963539
 2.05398  3.38753  2.51174  -0.996847     2.0  0.0  0.963539  0.963539
 2.01156  3.23139  2.37644  -0.714993  …  2.0  0.0  0.963539  0.963539
 2.01497  2.84762  2.48929  -1.15969      2.0  0.0  0.963539  0.963539
 2.05198  3.18417  2.31318  -0.711841     2.0  0.0  0.963539  0.963539
 2.49279  2.88198  2.45869  -1.05922      2.0  0.0  0.963539  0.963539
 2.01664  2.70077  2.65417  -0.983525     2.0  0.0  0.963539  0.963539
 2.33802  3.23896  2.5254   -1.05236   …  2.0  0.0  0.963539  0.963539
 1.73717  2.79669  2.31026  -0.831214     2.0  0.0  0.963539  0.963539
 2.16339  2.72781  2.37655  -0.925699     2.0  0.0  0.963539  0.963539
 ⋮                                     ⋱                      ⋮
 2.05396  3.58051  2.13081  -1.01232      2.0  0.0  0.963539  0.963539
 2.10231  3.03291  2.20647  -1.07199      2.0  0.0  0.963539  0.963539
 2.10231  3.03291  2.20647  -1.07199   …  1.0  0.0  0.963539  0.963539
 1.99741  3.06072  2.59206  -0.80927      2.0  0.0  0.963539  0.963539
 1.54651  3.0422   2.37786  -1.03687      2.0  0.0  0.963539  0.963539
 2.16404  3.1245   2.53922  -0.777052     2.0  0.0  0.963539  0.963539
 2.06172  2.946    2.29504  -1.04885      2.0  0.0  0.963539  0.963539
 1.92062  3.1204   2.52605  -0.775365  …  2.0  0.0  0.963539  0.963539
 2.14108  3.02875  2.51964  -1.02001      2.0  0.0  0.963539  0.963539
 1.94667  2.99099  2.28703  -0.833088     2.0  0.0  0.963539  0.963539
 2.08583  2.8571   2.53715  -0.91628      2.0  0.0  0.963539  0.963539
 2.08659  3.72413  2.3523   -0.85725      2.0  0.0  0.963539  0.963539
"c.value.data"
1000×16×2 Array{Float64, 3}:
[:, :, 1] =
 1.79171  2.9259   2.47615  -1.00088   …  2.0  0.0  0.943774  0.943774
 1.79171  2.9259   2.47615  -1.00088      2.0  0.0  0.943774  0.943774
 1.78365  3.20117  2.44422  -0.715697     2.0  0.0  0.943774  0.943774
 2.04188  2.71877  2.36068  -0.952244     2.0  0.0  0.943774  0.943774
 2.0008   3.34844  2.48828  -0.923913     2.0  0.0  0.943774  0.943774
 2.04423  2.64108  2.39721  -0.955044  …  2.0  0.0  0.943774  0.943774
 2.02032  3.60037  2.29273  -0.911933     2.0  0.0  0.943774  0.943774
 2.11052  2.62496  2.50005  -0.97364      2.0  0.0  0.943774  0.943774
 1.96377  3.46645  2.45473  -1.09452      2.0  0.0  0.943774  0.943774
 2.29862  3.644    2.55443  -0.889306     2.0  0.0  0.943774  0.943774
 1.76988  3.4902   2.38863  -0.820608  …  2.0  0.0  0.943774  0.943774
 1.94839  3.01624  2.55037  -0.704091     2.0  0.0  0.943774  0.943774
 2.07774  2.84975  2.38127  -0.993518     1.0  0.0  0.943774  0.943774
 ⋮                                     ⋱                      ⋮
 1.84218  3.05626  2.29374  -0.911148     2.0  0.0  0.943774  0.943774
 1.99603  3.36006  2.40845  -1.07951      2.0  0.0  0.943774  0.943774
 1.98988  2.53361  2.374    -0.787454  …  3.0  0.0  0.943774  0.943774
 1.98688  2.67816  2.35937  -1.17963      2.0  0.0  0.943774  0.943774
 1.74054  2.76433  2.4468   -0.91494      2.0  0.0  0.943774  0.943774
 1.6844   3.23351  2.39102  -0.881826     2.0  0.0  0.943774  0.943774
 1.70373  2.93048  2.36806  -1.00723      2.0  0.0  0.943774  0.943774
 1.95558  3.21013  2.29591  -0.785725  …  2.0  0.0  0.943774  0.943774
 2.25762  3.5584   2.14025  -0.684222     2.0  0.0  0.943774  0.943774
 2.02272  3.16815  2.28926  -0.661574     2.0  0.0  0.943774  0.943774
 1.98888  2.85139  2.28506  -1.10707      2.0  0.0  0.943774  0.943774
 2.22531  2.53992  2.42093  -0.805889     2.0  0.0  0.943774  0.943774

[:, :, 2] =
 1.72913  3.374    2.48625  -0.850005  …  2.0  0.0  0.963539  0.963539
 2.34423  2.83456  2.1521   -0.915034     2.0  0.0  0.963539  0.963539
 2.12179  3.28248  2.55795  -0.98139      2.0  0.0  0.963539  0.963539
 1.93434  2.76009  2.25582  -0.898026     2.0  0.0  0.963539  0.963539
 2.05398  3.38753  2.51174  -0.996847     2.0  0.0  0.963539  0.963539
 2.01156  3.23139  2.37644  -0.714993  …  2.0  0.0  0.963539  0.963539
 2.01497  2.84762  2.48929  -1.15969      2.0  0.0  0.963539  0.963539
 2.05198  3.18417  2.31318  -0.711841     2.0  0.0  0.963539  0.963539
 2.49279  2.88198  2.45869  -1.05922      2.0  0.0  0.963539  0.963539
 2.01664  2.70077  2.65417  -0.983525     2.0  0.0  0.963539  0.963539
 2.33802  3.23896  2.5254   -1.05236   …  2.0  0.0  0.963539  0.963539
 1.73717  2.79669  2.31026  -0.831214     2.0  0.0  0.963539  0.963539
 2.16339  2.72781  2.37655  -0.925699     2.0  0.0  0.963539  0.963539
 ⋮                                     ⋱                      ⋮
 2.05396  3.58051  2.13081  -1.01232      2.0  0.0  0.963539  0.963539
 2.10231  3.03291  2.20647  -1.07199      2.0  0.0  0.963539  0.963539
 2.10231  3.03291  2.20647  -1.07199   …  1.0  0.0  0.963539  0.963539
 1.99741  3.06072  2.59206  -0.80927      2.0  0.0  0.963539  0.963539
 1.54651  3.0422   2.37786  -1.03687      2.0  0.0  0.963539  0.963539
 2.16404  3.1245   2.53922  -0.777052     2.0  0.0  0.963539  0.963539
 2.06172  2.946    2.29504  -1.04885      2.0  0.0  0.963539  0.963539
 1.92062  3.1204   2.52605  -0.775365  …  2.0  0.0  0.963539  0.963539
 2.14108  3.02875  2.51964  -1.02001      2.0  0.0  0.963539  0.963539
 1.94667  2.99099  2.28703  -0.833088     2.0  0.0  0.963539  0.963539
 2.08583  2.8571   2.53715  -0.91628      2.0  0.0  0.963539  0.963539
 2.08659  3.72413  2.3523   -0.85725      2.0  0.0  0.963539  0.963539

Variable Types and Type Parameters¶

랜덤 변수의 벡터(또는 행렬)의 엘리먼트 유형은 사전분포의 eltype과 일치해야 합니다. 불연속 분포의 경우 <: Integer, 연속 분포의 경우 <: AbstractFloat여야 합니다.또한 해밀턴 샘플러를 사용하여 연속형 확률 변수를 샘플링하려면 벡터의 요소 유형이 1. 실수(Real)여야 합니다. 실수(Real)의 하위 유형인 특정 숫자 유형을 사용하는 모델을 통해 자동미분(auto-diffrentiation) 이 가능하거나 2. type parameter 구문을 사용하여 모델 헤더에 정의된 일부 type parameter T(예: function gdemo(x, ::Type{T} = Float64) where {T}). 마찬가지로 파티클 샘플러를 사용할 때 사용되는 Julia 변수는 1. TArray 이거나 2. type parameter 구문을 사용하여 모델 헤더에 정의된 type parameter T의 인스턴스(예: function gdemo(x, ::Type{T} = Vector{Float64}) where {T}).

Querying Probabilities from Model or Chain¶

  • `prob` macro 사용시 model의 파라미터는 model의 instance가 되어야 한다.
  • 예) model_instance = gdemo(0,0) gdemo의 입력값은 아무거나
  • likelihood : $p(x=2.5|s=1.2)$ = prob"x=2.5 , s=1.2| model = m"
  • posterior : $p(s=1.2|x=2.5)$ = prob"x=2.5 | model=m, s=1.2"
  • prior: $p(s=1.2)$ = prob"s=1.2 | model = m, x = nothing"
In [22]:
let
    @model function gdemo(x)
        s ~ Normal(0, 1)
        x ~ Normal(s, 1)
    end
    
    # model instance를 만들기 위해서 입력값을 아무거나 넣어도 되지만
    # nothing 으로 설정
    m = gdemo(nothing)
    
    
    # 사전확률밀도 계산
    prior_01 = prob"s=1.2 | model = m, x = nothing"
    # 직접계산
    prior_02 = pdf(Normal(0,1), 1.2)
    
    println("prior_01 : ", prior_01)
    println("prior_02 : ", prior_02)
    
    # 우도 계산
    likelihood_01 = prob"x=2.5 |model = m, s=1.2"
    # 직접계산
    # s = 1.2, x=2.5 
    likelihood_02 = pdf(Normal(1.2,1),2.5)
    
    println("likelihood_01 : ", likelihood_01)
    println("likelihood_02 : ", likelihood_02)
    
    # 사후확률밀도 계산
    # 사후확률밀도 = 우도 * 사전확률밀도
    posterior_01 = prob"s=1.2, x=2.5 | model=m"
    # 직접계산
    posterior_02 =  pdf(Normal(1.2,1),2.5)*pdf(Normal(0,1),1.2)
    
    println("posterior_01: ", posterior_01)
    println("posterior_02: ", posterior_02)
end

prior_01 : 0.19418605498321292
prior_02 : 0.19418605498321295
likelihood_01 : 0.17136859204780733
likelihood_02 : 0.17136859204780736
posterior_01: 0.0332773908377913
posterior_02: 0.03327739083779131
In [23]:
let
    using Plots, KernelDensity   
    plotlyjs()   

    # x,y : i.i.d
    # m,s²: i.i.d
    
    
    @model function gdemo(x,y)
        s²~ InverseGamma(2, 3)
        m ~ Normal(0, √s²)
        x ~ Normal(m, √s²)
        y ~ Normal(m, √s²)
    end    
    
    model_instance= gdemo(nothing,nothing)
    
    pdf_s2(v) = pdf(InverseGamma(2,3),v)

    #--------------------------------
    # Likelihood : x,y
    #--------------------------------
    likelihood_01 = prob"x = 1.0, y = 3.0 | s² = 1.5, m = 2.0, model = model_instance"     
    # p(x,y|m,s²) = p(x|m,s²)*p(y|m,s²)
    likelihood_02 = pdf(Normal(2.0,√1.5),1.0) * pdf(Normal(2.0,√1.5),3.0)
    
    println("likehood_01 : ",likelihood_01)
    println("likehood_02 : ",likelihood_02)
    
    #--------------------------------
    # Joint probability : m,s²
    #--------------------------------
    joint_prob1_01 = prob"s² = 1.5, m = 2.0 | x = nothing, y = nothing, model = model_instance"         
    # p(s²,m) = p(s²) *  p(m)
    joint_prob1_02 = pdf_s2(1.5)*pdf(Normal(0,√1.5),2.0)
    
    println("joint_prob1_01 of m, s² : ",joint_prob1_01)
    println("joint_prob1_02 of m, s² : ",joint_prob1_02)
    
    #--------------------------------
    # Joint probability : m,s²,x
    #--------------------------------
    joint_prob2_01 = prob"s² = 1.5, m = 2.0, x = 1.0 | y = nothing, model = model_instance"         
    # p(s²,m,x) = p(s²) *  p(m) * p(x)
    joint_prob2_02 = pdf_s2(1.5)*pdf(Normal(0,√1.5),2.0)*pdf(Normal(2.0,√1.5),1.0)
    
    println("joint_prob2_01 of m, s²,x : ",joint_prob2_01)
    println("joint_prob2_02 of m, s²,x : ",joint_prob2_02)

    #--------------------------------
    # Joint probability : m,s²,x, y
    #--------------------------------
    joint_prob3_01 = prob"s² = 1.5, m = 2.0, x = 1.0, y = 3.0 | model = model_instance"         
    # p(s²,m,x) = p(s²) *  p(m) * p(x)
    joint_prob3_02 = pdf_s2(1.5)*pdf(Normal(0,√1.5),2.0)*
                     pdf(Normal(2.0,√1.5),1.0)*pdf(Normal(2.0,√1.5),3.0)
    
    println("joint_prob3_01 of m, s²,x, y : ",joint_prob3_01)
    println("joint_prob3_02 of m, s²,x, y : ",joint_prob3_02)    
    
    #--------------------------------
    # MCMC Sampling 이후
    #--------------------------------   
    chain = sample(gdemo(1.0, 3.0), NUTS(), 1000)
    plot(chain) |> display
    autocorplot(chain) |> display
    
    # DataFrame으로 chian 데이터 읽어 오기
    df = DataFrame(chain)
    df[1:5,:] |> display
    
    #--------------------------------
    # Likelihood : x,y
    # MCMC에서 샘플링한 1000개의 m, s²에 대해 
    #  x,y의 likelihood 계산
    #  p(x,y | m,s²)
    #--------------------------------
    # 경고 숨기기
    likelihood_mcmc = @suppress begin
        likelihood_mcmc = prob"x = 1.0, y = 3.0 | chain=chain, model = model_instance"
        likelihood_mcmc
    end    

    println("likelihood_mcmc size : ",size(likelihood_mcmc))
    println("size(chain[:s²]) : ",size(chain[:s²]))
    println("size(chain[:m]) : ",size(chain[:m]))

    #--------------------------------
    # Sampleing된 m, s²데이터로 부터 
    # density를 구하고 argmax를 구한다
    # KDE 방법으로 데이터 밀도 추정
    #--------------------------------
    s²_density = kde(df.s²)
    s²_max_idx = argmax(s²_density.density)    
    s²_map = s²_density.x[s²_max_idx]
    
    m_density = kde(df.m)
    m_max_idx = argmax(m_density.density)    
    m_map = m_density.x[m_max_idx]
    
    #plot(m_density.x,m_density.density) |> display
    
    #-------------------------------
    # Visualize s²,m with KDE
    #-------------------------------
    plot(s²_density.x, s²_density.density, label="", lw=2)
    vline!([s²_map], color=:red, label="Max Density at s²= $(round(s²_map,digits=4))",
            linestyle=:dash,lw=2)
    xlabel!("s² value")
    ylabel!("Density")
    p1 = title!("Kernel Density Estimation of s²")
    
    plot(m_density.x, m_density.density, label="", lw=2)
    vline!([m_map], color=:red, label="Max Density at m = $(round(m_map,digits=4))", lw=2)
    xlabel!("m value")
    ylabel!("Density")
    p2 = title!("Kernel Density Estimation of m")

    plot(p1,p2, layout= (1,2), legend=:bottomright, titlefontsize=12) |> display
    
    
    println("Max Posterior of s²,m : $(s²_map), $(m_map)")
    
    #--------------------------------
    # Visualize
    #--------------------------------    
    # (1000,1) => (1000,) 으로 변환
    l_mcmc = reshape(likelihood_mcmc,(:,))
    idx = argmax(l_mcmc)
    println("max likelihood index : ",idx)
    df[idx:idx,:] |> display
    println("max likelihood : p(x=1.0, y=3.0 | s²= $(df.s²[idx]),m=$(df.m[idx])) = $(l_mcmc[idx]) ") 

    sp1 = scatter3d(df.s²,df.m, l_mcmc, markersize=0.3,legend=false,
        xlabel="s²", ylabel="m", zlabel="p(x=1.0,y=3.0 | s²,m)", size=(1000,800))
    sp2 = scatter3d(sp1,[df.s²[idx]],[df.m[idx]], [l_mcmc[idx]], markersize=1.0,color=:red) 
    x,y,z = df.s²[idx],df.m[idx],l_mcmc[idx]
    x_line = [x, x, x]
    y_line = [0, y, y]
    z_line = [0, 0, z]
    plot3d!(sp2,x_line, y_line,z_line, mode="lines", color=:red, linewidth=5)
end

The WebIO Jupyter extension was not detected. See the WebIO Jupyter integration documentation for more information. 

likehood_01 : 0.05447524824135803
likehood_02 : 0.054475248241358035
joint_prob1_01 of m, s² : 0.030987382682801004
joint_prob1_02 of m, s² : 0.03098738268280102
joint_prob2_01 of m, s²,x : 0.007232434422796936
joint_prob2_02 of m, s²,x : 0.007232434422796939
joint_prob3_01 of m, s²,x, y : 0.0016880453639955442
joint_prob3_02 of m, s²,x, y : 0.0016880453639955449

┌ Info: Found initial step size
└   ϵ = 0.8
5×16 DataFrame
Rowiterationchains²mlpn_stepsis_acceptacceptance_ratelog_densityhamiltonian_energyhamiltonian_energy_errormax_hamiltonian_energy_errortree_depthnumerical_errorstep_sizenom_step_size
Int64Int64Float64Float64Float64Float64Float64Float64Float64Float64Float64Float64Float64Float64Float64Float64
150113.202740.451147-6.663353.01.00.995849-6.663357.38575-0.0443729-0.2972672.00.00.8567250.856725
250211.695260.855131-5.755383.01.00.996916-5.755386.71051-0.128104-0.1281042.00.00.8567250.856725
350315.2232-0.481347-8.312277.01.00.790509-8.312278.43730.3176790.9726322.00.00.8567250.856725
450415.591030.646362-7.664183.01.00.938712-7.664188.994530.0632611-0.3025262.00.00.8567250.856725
550514.17639-0.646191-8.247063.01.01.0-8.247068.81888-0.071352-0.4590792.00.00.8567250.856725
likelihood_mcmc size : (1000, 1)
size(chain[:s²]) : (1000, 1)
size(chain[:m]) : (1000, 1)

Max Posterior of s²,m : 1.4815249676020907, 1.2829839545006743
max likelihood index : 554
1×16 DataFrame
Rowiterationchains²mlpn_stepsis_acceptacceptance_ratelog_densityhamiltonian_energyhamiltonian_energy_errormax_hamiltonian_energy_errortree_depthnumerical_errorstep_sizenom_step_size
Int64Int64Float64Float64Float64Float64Float64Float64Float64Float64Float64Float64Float64Float64Float64Float64
1105410.9285372.01965-6.804813.01.00.877812-6.804817.406480.37980.37982.00.00.8567250.856725
max likelihood : p(x=1.0, y=3.0 | s²= 0.9285366482105682,m=2.019651684307102) = 0.05836079317783932
Out[23]:

Maximum likelihood and maximum a posterior estimates¶

Turing은 두 가지 모드 추정 기법, 즉 최대우도추정(MLE-maximum likelihood estimation)과 최대사후추정(MAP-maximum a posterior)을 지원합니다. 최적화는 Optim.jl 패키지에 의해 수행됩니다.

  • MLE의 결과를 위의 sample에서 구한 값과 비교 해볼것
  • MAP의 결고를 위의 sample에서 구한 값과 비교 해볼것

일부 메서드는 충분한 반복이 허용되지 않거나 함수 호출 사이에 대상 함수가 위로 이동하여 모드를 계산하는 데 문제가 있을 수 있습니다. Optim.converge를 실행하여 Optim이 수렴에 실패하면 Turing에서 경고를 표시합니다. 이에 대한 일반적인 해결책은 반복을 더 추가하거나 함수 반복 사이에 옵티마이저가 증가하도록 허용하는 것입니다:

In [24]:
(mle_estimate, map_estimate,model) = let
    using Optim
    
    @model function gdemo(x)
        s²~InverseGamma(2,3)
        m ~ Normal(0, √s²)
        
        for i in eachindex(x)
            x[i] ~ Normal(m, √s²)
        end
    end
    
    data = [1.0, 3.0]
    
    model = gdemo(data)
    
    #---------------------------------
    # MLE, MAP Optimization
    # default optimizer : LBFGS optimizer
    #---------------------------------
    # Generate a MLE estimate.
    mle_estimate = optimize(model, MLE())
    
    # Generate a MAP estimate.
    map_estimate = optimize(model, MAP())
        
    fieldnames(typeof(map_estimate)) |> display
    
    println("MLE of x=1.0, y=3.0 : ")
    mle_estimate |> display
    
    println("MAP of x=1.0, y=3.0 : ")
    map_estimate |> display
    #---------------------------------
    # Use NelderMead
    #---------------------------------
    # Increase the iterations and allow function eval to increase between calls.
    mle_estimate2 = optimize(model, MLE(), NelderMead())
    println("MLE-NelderMead of x=1.0, y=3.0 : ")
    mle_estimate2 |> display         
    
    (mle_estimate, map_estimate,model)
end;

(:values, :optim_result, :lp, :f)
MLE of x=1.0, y=3.0 : 

ModeResult with maximized lp of -2.84
[1.000000000710164, 2.0000000009077112]
MAP of x=1.0, y=3.0 : 

ModeResult with maximized lp of -5.82
[1.1851851851851711, 1.3333333333333324]
MLE-NelderMead of x=1.0, y=3.0 : 

ModeResult with maximized lp of -2.84
[0.999937799356961, 1.999951237477533]

Analyzing your mode estimate¶

Turing은 모드 추정 결과를 분석하는 데 사용할 수 있는 StatsBase의 여러 메서드를 확장합니다. 구현된 메서드에는 vcov, informationmatrix, coeftable, params, coef 등이 있습니다.

표준 오차는 피셔 정보 행렬(로그 확률 또는 로그 조인트의 역 헤시안)에서 계산되며, 빈도 통계학자에게는 t-통계가 친숙할 것입니다. 경고 - 이러한 방식으로 계산된 표준 오차는 MAP 추정치에 항상 적합하지 않을 수 있으므로 해석에 주의하세요.

For example, let's examine our ML estimate from above using coeftable:

In [25]:
let (mle_estimate, map_estimate) = (mle_estimate, map_estimate)
    using StatsBase
    
    # Print out the coefficient table
    coeftable(mle_estimate) |> DataFrame |> display
    coeftable(map_estimate) |> DataFrame |> display
    
end
2×7 DataFrame
RowNameCoef.Std. ErrorzPr(>|z|)Lower 95%Upper 95%
StringFloat64Float64Float64Float64Float64Float64
1s²1.01.01.00.317311-0.9599642.95996
2m2.00.7071072.828430.004677730.6140963.3859
2×7 DataFrame
RowNameCoef.Std. ErrorzPr(>|z|)Lower 95%Upper 95%
StringFloat64Float64Float64Float64Float64Float64
1s²1.185190.5587022.121320.03389490.09015012.28022
2m1.333330.6285392.121320.03389490.1014192.56525

Sampling with the MAP/MLE as initial states¶

파라미터 값의 벡터를 추출하고 키워드 init_params를 사용하여 샘플 함수에 제공함으로써 MLE/MAP 추정치에서 체인 샘플링을 시작할 수 있습니다. 예를 들어, 다음은 MAP 추정치를 시작점으로 사용하여 전체 사후분포(posterior)에서 샘플링하는 방법입니다:

In [26]:
let model = model
    # G
    map_estimate =optimize(model, MAP())
    
    coeftable(map_estimate) |> DataFrame |> display
    
    # Sample with the MAP estimate as the starting point.
    chain = sample(model, NUTS(), 1_000, init_params = map_estimate.values.array)
    #chain |> DataFrame |> display
    plot(chain,titlefontsize=10) |> display
end
2×7 DataFrame
RowNameCoef.Std. ErrorzPr(>|z|)Lower 95%Upper 95%
StringFloat64Float64Float64Float64Float64Float64
1s²1.185190.5587022.121320.03389490.09015012.28022
2m1.333330.6285392.121320.03389490.1014192.56525
┌ Info: Found initial step size
└   ϵ = 1.6
Sampling: 100%|█████████████████████████████████████████| Time: 0:00:00

Beyond the Basic¶

Compositional Sampling Using Gibbs¶

Turing.jl은 서로 다른 샘플러를 결합할 수 있는 Gibbs 인터페이스를 제공합니다. 예를 들어, 아래와 같이 HMC 샘플러와 PG 샘플러를 결합하여 단일 모델에서 다양한 파라미터에 대한 추론을 실행할 수 있습니다.

깁스 샘플러를 사용하여 다양한 변수 공간에 대해 고유한 자동미분(Automatic Differentiation) 백엔드를 지정할 수 있습니다. 자세한 내용은 자동미분(Automatic Differentiation) 문서를 참조하세요.

Turing.jl의 컴포지션 샘플링에 대한 자세한 내용은 해당 논문을 참조하세요.

In [27]:
let
    @model function simple_choice(xs)
        p ~ Beta(2, 2)
        z ~ Bernoulli(p)
        for i in 1:length(xs)
            if z == 1
                xs[i] ~ Normal(0, 1)
            else
                xs[i] ~ Normal(2, 1)
            end
        end
    end

    simple_choice_f = simple_choice([1.5, 2.0, 0.3])

    chn = sample(simple_choice_f, Gibbs(HMC(0.2, 3, :p), PG(20, :z)), 1000)
    plot(chn) |> display
    chn |> display
end

Sampling: 100%|█████████████████████████████████████████| Time: 0:00:16

Chains MCMC chain (1000×3×1 Array{Float64, 3}):

Iterations        = 1:1:1000
Number of chains  = 1
Samples per chain = 1000
Wall duration     = 27.65 seconds
Compute duration  = 27.65 seconds
parameters        = p, z
internals         = lp

Summary Statistics
  parameters      mean       std      mcse   ess_bulk   ess_tail      rhat   e ⋯
      Symbol   Float64   Float64   Float64    Float64    Float64   Float64     ⋯

           p    0.3929    0.2087    0.0329    42.9175   176.4490    1.0389     ⋯
           z    0.1270    0.3331    0.0162   425.3284        NaN    0.9990     ⋯
                                                                1 column omitted

Quantiles
  parameters      2.5%     25.0%     50.0%     75.0%     97.5% 
      Symbol   Float64   Float64   Float64   Float64   Float64 

           p    0.0694    0.2252    0.3734    0.5460    0.8089
           z    0.0000    0.0000    0.0000    0.0000    1.0000

Working with filldist and arraydist¶

Turing은 production 분포를 구성하기 위한 단순화된 인터페이스로 filldist(dist::Distribution, n::Int) 및 arraydist(dists::AbstractVector{<:Distribution})를 제공합니다(예: 동일한 구조를 공유하지만 그룹마다 다른 변수 집합을 모델링하기 위해).

a ~ filldist(Exponential(), k) 는 아래와 동일

a[1] ~ Exponential()

...

a[k] ~ Exponential()

Constructing product disttributions with filldist¶

filldist 함수는 동일한 유형 및 매개변수화 분포에 대한 분포곱하기을 구성하는 일반 인터페이스를 제공합니다. Distributions.jl(Product)에서 제공하는 Production Distribution 인터페이스와 달리 filldist는 단변량 또는 다변량 분포에 대한 분포 곱하기를 지원한다는 점에 유의하세요.

**코드설명(by Chat-GPT4)**

아래 코드는 주어진 x 데이터와 그룹 정보 g를 사용하여 간단한 계층적 모델을 정의하는 Julia 코드입니다. 모델의 주요 개념과 구성 요소를 설명하겠습니다:

  1. 함수 정의: @model 매크로는 Turing.jl에서 확률 모델을 정의할 때 사용합니다. demo라는 함수가 정의되었으며, 이 함수는 x (관측 데이터)와 g (그룹 또는 카테고리 정보) 두 가지 인수를 받습니다.

  2. 그룹의 수 계산: k = length(unique(g)) - g 내의 고유한 그룹 또는 카테고리의 수를 계산합니다. 예를 들어, g = [1, 1, 2, 3, 3]이라면 k는 3이 됩니다.

  3. 계층적 사전 분포: a ~ filldist(Exponential(), k) - k 개의 그룹 각각에 대해 지수 분포 (Exponential())로부터 뽑힌 값을 갖는 a 벡터를 정의합니다. 즉, 각 그룹에 대한 파라미터 a를 독립적으로 지수 분포에서 샘플링합니다.

  4. 모델링 mu: mu = a[g] - 데이터 포인트 각각에 대한 평균 값을 나타내는 mu를 계산합니다. 여기서 g의 각 원소는 해당 데이터 포인트가 속하는 그룹을 나타내므로, mu는 해당 그룹의 a 값으로 설정됩니다.

  5. 관측 데이터에 대한 가능도: x .~ Normal.(mu) - 관측 데이터 x 각각은 해당 데이터 포인트에 대한 mu 값으로 평균을 갖는 정규 분포에서 독립적으로 샘플링되었다고 가정합니다.

이 모델은 그룹화 된 데이터에 대한 계층적 구조를 가진 확률 모델을 나타냅니다. 각 그룹은 독립적인 지수 분포에서 파라미터를 가지며, 각 데이터 포인트는 해당 그룹의 파라미터를 평균으로 하는 정규 분포에서 샘플링된 것으로 가정합니다.

In [28]:
let
    @model function demo(x, g)
        k = length(unique(g))
        a ~ filldist(Exponential(), k) # = Product(fill(Exponential(), k))
        mu = a[g]
        x .~ Normal.(mu)
    end
    x = [rand(Normal(0,1),30);] #rand(Normal(3,1),3);rand(Normal(5,1),3);]
    g = rand(1:3,length(x))
    model = demo(x,g)
    samples = sample(model,NUTS(),1000)
    plot(samples, titlefontsize=10)
end

┌ Info: Found initial step size
└   ϵ = 0.05
Sampling: 100%|█████████████████████████████████████████| Time: 0:00:00
Out[28]:

Constructing product distributions with arraydist¶

arraydist 함수는 다양한 유형 및 매개변수화 분포에 대한 production 분포를 구성하는 일반적인 인터페이스를 제공합니다. Distributions.jl(Product)에서 제공하는 Production 분포 인터페이스와 달리 arraydist는 단변량 또는 다변량 분포에 대한 Production 분포를 지원한다는 점에 유의하세요.

In [29]:
let
    @model function demo(x, g)
        k = length(unique(g))
        a ~ arraydist([Exponential(i) for i in 1:k])
        mu = a[g]
        x .~ Normal.(mu)
    end
    x = [rand(Normal(0,1),30);] #rand(Normal(3,1),3);rand(Normal(5,1),3);]
    g = rand(1:3,length(x))
    model = demo(x,g)
    samples = sample(model,NUTS(),1000)
    plot(samples, titlefontsize=10)
end

┌ Info: Found initial step size
└   ϵ = 0.2
Sampling: 100%|█████████████████████████████████████████| Time: 0:00:00
Out[29]:

계층적 확률 모델 이란?¶

계층적 확률 모델은 데이터의 계층적 구조를 반영하기 위해 설계된 모델입니다. 이러한 모델은 다수의 관측 단위와 그룹 간 및 그룹 내 변동을 설명하는 변수들을 포함합니다. 계층적 모델은 특히 그룹화된 데이터에서 유용하며, 각 그룹이 고유한 패턴을 가질 수 있지만 전체적으로 일관된 추세를 공유하는 경우에 많이 사용됩니다.

기본 개념¶

  1. 계층적 모델링: 데이터에 내재된 계층 구조를 반영하여 파라미터를 추정합니다. 각 계층은 그룹별 차이와 전체적인 차이를 모두 포함합니다.

  2. 분산 구조: 그룹 내와 그룹 간 변동성을 모두 모델링하여, 각 그룹의 특성과 전체 특성을 모두 반영합니다.

예제: 학교별 학생 성적¶

학생들의 시험 점수를 모델링하려고 할 때, 각 학생은 특정 학교에 속합니다. 학생들의 점수는 그들의 개인적인 능력 뿐만 아니라 학교의 특성에 따라 다를 수 있습니다.

  1. 첫 번째 계층: 각 학생의 개별 성적을 모델링합니다.
  2. 두 번째 계층: 학교별 평균 점수와 변동성을 모델링합니다.

모델은 일반적으로 다음과 같이 나타낼 수 있습니다:

  • $y_{ij}$ = i번째 학교의 j번째 학생의 점수
  • $\alpha$ = 전체 학교의 평균 점수
  • $\sigma$ = 학교간 평균점수간의 표준편차
  • $\mu_i$ = i번째 학교의 평균 점수, $\mu_i \sim \mathcal{N}(\alpha,\sigma^2)$
  • $\epsilon_{ij}$ = i번째 학교의 평균점수와 j번째 학생의 점수간의 표준편차

모델은 다음과 같이 표현됩니다:
$y_{ij} = \mu_i + \epsilon_{ij}$

여기서, ( $\mu_i$ )는 학교 효과를 나타내며, 이 효과는 전체 학교에 대한 평균 효과 ( $\alpha$ )와 함께 모델에 포함됩니다.

이러한 계층적 구조를 사용하면 각 학교의 특성과 전체적인 특성을 동시에 추정할 수 있습니다. 그 결과, 특정 학교의 데이터가 부족한 경우에도 전체적인 추세를 기반으로 그 학교에 대한 추정을 더 정확하게 할 수 있습니다.

학교별 학생 성적을 모델링하는 간단한 Turing 예제 (by Chat-GPT 4)¶

먼저, 간단한 예제 데이터를 생성합니다:

5개의 학교에서 각각 30명의 학생들의 성적을 생성합니다. 각 학교마다 고유한 평균 점수를 가정하며, 학생들의 점수는 학교의 평균을 중심으로 랜덤하게 분포됩니다.

이 모델에서:

  • α는 전체 평균 점수를 나타냅니다.
  • σ는 학교 간 점수의 분산을 나타냅니다.
  • μ는 각 학교별 평균 점수를 나타냅니다.
  • 각 학생의 점수는 그 학교의 평균을 중심으로 정규분포를 따릅니다.

NUTS() sampler를 사용하여 모델을 샘플링하고, 결과를 chain 변수에 저장합니다.

아래 코드에서의 확률 모델을 수식으로 나타내면 다음과 같습니다:¶

사전분포 (Prior Distributions):¶

  1. 전체 학교 평균에 대한 사전분포:

$ \alpha \sim \mathcal{N}(0, 10^2) $

  1. 학교 간 분산에 대한 사전분포:

$ \sigma \sim \text{Truncated}(\text{Cauchy}(0, 10), 0, \infty) $

  1. 학교별 평균에 대한 사전분포:

$ \mu_j \sim \mathcal{N}(\alpha, \sigma^2) \quad \text{for } j = 1, \ldots, n_{\text{schools}} $

우도 (Likelihood):¶

학생별 점수에 대한 우도: $ \epsilon = 5\\ \text{grades}_i \sim \mathcal{N}(\mu_{\text{school_ids}[i]}, \epsilon^2) \quad \text{for } i = 1, \ldots, \text{length(grades)} \\ $

사후분포 (Posterior Distributions):¶

사후분포는 위의 사전분포와 우도를 조합하여 계산됩니다. 수식으로 간단하게 표현하면:

$ p(\alpha, \sigma, \mu | \text{grades}) \propto p(\alpha) p(\sigma) \prod_{j=1}^{n_{\text{schools}}} p(\mu_j | \alpha, \sigma) \prod_{i=1}^{\text{length(grades)}} p(\text{grades}_i | \mu_{\text{school_ids}[i]}) $

여기서, $p(\alpha)$ , $p(\sigma)$ , $p(\mu_j | \alpha, \sigma)$ , 그리고 $p(\text{grades}_i | \mu_{\text{school_ids}[i]})$ 는 각각 위에서 설명한 사전분포와 우도를 나타냅니다.

filldist(Normal(α, σ), 3)¶

μ ~ filldist(Normal(α, σ), 3) 는 아래와 같은 내용임

μ[1] ~ Normal(α, σ), μ[2] ~ Normal(α, σ), μ[2] ~ Normal(α, σ)

결과 해석¶

전체학교 평균 = $\alpha + \sigma$

In [30]:
st_chain = let
    using Turing, Distributions, Random

    # 예제 데이터 생성
    Random.seed!(1234) # 재현성을 위한 시드 설정

    n_schools = 5
    students_per_school = 30

    # 학교별 평균 점수
    school_avg = [85, 90, 75, 80, 70]

    # 각 학교의 학생들의 점수 생성
    grades = []
    school_ids = []
    for i in 1:n_schools
        school_grades = rand(Normal(school_avg[i], 10), students_per_school)
        push!(grades, school_grades)
        push!(school_ids, fill(i, students_per_school))
    end
    grades = vcat(grades...)
    school_ids = vcat(school_ids...)

    # 모델 정의
    @model function hierarchical_model(grades, school_ids, n_schools)
        α ~ Normal(0, 10) # 전체 학교 평균에 대한 사전분포 (0점이라고 가정 - 무리한 가정)
        σ ~ Truncated(Cauchy(0, 10), 0, Inf) # 학교 간 분산
        μ ~ filldist(Normal(α, σ), n_schools) # 학교별 평균
        for i in 1:length(grades)
            grades[i] ~ Normal(μ[school_ids[i]], 5) # 학생별 점수
        end
    end

    # 샘플링
    chain = sample(hierarchical_model(grades, school_ids, n_schools), NUTS(), 1000)
    plot(Truncated(Cauchy(0, 10), 0, 200), title="σ ~ Truncated(Cauchy(0, 10), 0, Inf)",
        titlefontsize=10, legend=false, size=(800,200)) |> display
    plot(chain, titlefontsize=8) |> display
    chain
end

┌ Info: Found initial step size
└   ϵ = 0.8
Sampling: 100%|█████████████████████████████████████████| Time: 0:00:00
Out[30]:
Chains MCMC chain (1000×19×1 Array{Float64, 3}):

Iterations        = 501:1:1500
Number of chains  = 1
Samples per chain = 1000
Wall duration     = 7.0 seconds
Compute duration  = 7.0 seconds
parameters        = α, σ, μ[1], μ[2], μ[3], μ[4], μ[5]
internals         = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size

Summary Statistics
  parameters      mean       std      mcse    ess_bulk   ess_tail      rhat    ⋯
      Symbol   Float64   Float64   Float64     Float64    Float64   Float64    ⋯

           α    8.9753   10.7545    0.4211    656.0845   555.2656    1.0051    ⋯
           σ   74.6378   28.7578    1.0974    561.8351   662.6308    1.0032    ⋯
        μ[1]   85.7657    0.9061    0.0240   1412.4900   444.3103    1.0038    ⋯
        μ[2]   89.6041    0.9454    0.0257   1360.8753   693.5409    1.0006    ⋯
        μ[3]   71.1830    0.9059    0.0305    886.0895   592.2101    1.0001    ⋯
        μ[4]   82.0604    0.9339    0.0258   1379.9769   611.3171    1.0069    ⋯
        μ[5]   70.4427    0.9181    0.0318    831.5316   679.7763    1.0034    ⋯
                                                                1 column omitted

Quantiles
  parameters       2.5%     25.0%     50.0%     75.0%      97.5% 
      Symbol    Float64   Float64   Float64   Float64    Float64 

           α   -11.4631    1.3864    8.4896   16.4780    30.2636
           σ    37.6910   54.8507   68.6761   87.8000   146.9226
        μ[1]    84.0517   85.1536   85.7371   86.3947    87.5047
        μ[2]    87.8123   88.9465   89.5805   90.2405    91.5007
        μ[3]    69.3960   70.5671   71.2126   71.7931    72.9761
        μ[4]    80.1276   81.4287   82.0722   82.7079    83.7709
        μ[5]    68.6682   69.8429   70.4582   71.0207    72.2448

Working with MCMCChain.jl¶

Turing.jl은 MCMCChains.Chain을 사용하여 샘플을 래핑하므로 MCMCChains.Chain에서 작동하는 모든 함수를 Turing.jl에서 재사용할 수 있습니다. 대표적인 두 가지 함수는 MCMCChains.describe와 MCMCChains.plot이며, 획득한 체인 chn에 대해 다음과 같이 사용할 수 있습니다. MCMCChains에 대한 자세한 내용은 GitHub 리포지토리를 참조하세요.

In [31]:
let
    describe(st_chain) |> display
    describe(st_chain["α"]) |> display
    describe(st_chain["σ"]) |> display
    describe(st_chain["μ[1]"]) |> display
end;

2-element Vector{ChainDataFrame}:
 Summary Statistics (7 x 8)
 Quantiles (7 x 6)
Summary Stats:
Length:         1000
Missing Count:  0
Mean:           8.975331
Minimum:        -22.544918
1st Quartile:   1.386443
Median:         8.489586
3rd Quartile:   16.478041
Maximum:        41.300121
Type:           Float64

nothing
Summary Stats:
Length:         1000
Missing Count:  0
Mean:           74.637763
Minimum:        28.109327
1st Quartile:   54.850704
Median:         68.676125
3rd Quartile:   87.799981
Maximum:        257.789839
Type:           Float64

nothing
Summary Stats:
Length:         1000
Missing Count:  0
Mean:           85.765668
Minimum:        82.353224
1st Quartile:   85.153615
Median:         85.737066
3rd Quartile:   86.394738
Maximum:        89.386477
Type:           Float64

nothing

Working with Libtask.jl¶

Libtask.jl 라이브러리는 Turing의 particle(입자) 기반 샘플러에서 사용하기에 안전한 write-on-copy 데이터 구조를 제공합니다. 특히 TArray를 사용해야 하는 경우가 많습니다. 다음 샘플러 유형은 분포를 저장하기 위해 TArray를 사용해야 합니다:

  • IPMCMC
  • IS
  • PG
  • PMMH
  • SMC

particle(입자) 기반 샘플러를 사용할 때 분포 배열을 저장할 때 TArray를 사용하지 않으면 오류가 발생할 수 있습니다.

In [32]:
let
    @model function BayesHmm(y,N, K)
        # Declare a TArray with a length of N.
        s = tzeros(Int, N)
        m = Vector{Float64}(undef, K)
        T = Vector{Vector{Float64}}(undef, K)
        for i = 1:K
            T[i] ~ Dirichlet(ones(K)/K)
            m[i] ~ Normal(i, 0.01)
        end
        
        # Draw from a distribution for each element in s.
        s[1] ~ Categorical(K)
        for i=2:N
            s[i] ~ Categorical(vec(T[s[i-1]]))
            y[i] ~ Normal(m[s[i]], 0.1)
        end
        return (s,m)
    end
    N = 3
    K = 2
    y = rand(Normal(3,1),N)
    model = BayesHmm(y,N,K)
    samples = sample(model, PG(K), 1000)
    plot(samples, titlefontsize=9) |> display
    samples
end

Sampling: 100%|█████████████████████████████████████████| Time: 0:00:03
Out[32]:
Chains MCMC chain (1000×11×1 Array{Float64, 3}):

Log evidence      = -293.45395345351056
Iterations        = 1:1:1000
Number of chains  = 1
Samples per chain = 1000
Wall duration     = 9.25 seconds
Compute duration  = 9.25 seconds
parameters        = T[1][1], T[1][2], T[2][1], T[2][2], m[1], m[2], s[1], s[2], s[3]
internals         = lp, logevidence

Summary Statistics
  parameters      mean       std      mcse   ess_bulk   ess_tail      rhat   e ⋯
      Symbol   Float64   Float64   Float64    Float64    Float64   Float64     ⋯

     T[1][1]    0.5438    0.4116    0.1093    23.5815    48.3905    1.0505     ⋯
     T[1][2]    0.4562    0.4116    0.1093    23.5815     8.5954    1.0505     ⋯
     T[2][1]    0.2529    0.2524    0.0675    20.0704    25.1402    1.0441     ⋯
     T[2][2]    0.7471    0.2524    0.0675    20.0704     8.9894    1.0441     ⋯
        m[1]    1.0042    0.0073    0.0015    22.6124    50.2898    1.0651     ⋯
        m[2]    2.0164    0.0052    0.0013    17.0329    29.6524    1.0903     ⋯
        s[1]    1.7070    0.4554    0.0792    33.0317        NaN    0.9991     ⋯
        s[2]    2.0000    0.0000       NaN        NaN        NaN       NaN     ⋯
        s[3]    1.9920    0.0891    0.0101    77.5916        NaN    1.0071     ⋯
                                                                1 column omitted

Quantiles
  parameters      2.5%     25.0%     50.0%     75.0%     97.5% 
      Symbol   Float64   Float64   Float64   Float64   Float64 

     T[1][1]    0.0005    0.1448    0.6160    0.9715    0.9904
     T[1][2]    0.0096    0.0285    0.3840    0.8552    0.9995
     T[2][1]    0.0004    0.0203    0.1608    0.4824    0.6642
     T[2][2]    0.3358    0.5176    0.8392    0.9797    0.9996
        m[1]    0.9867    0.9986    1.0066    1.0083    1.0137
        m[2]    2.0055    2.0136    2.0165    2.0183    2.0252
        s[1]    1.0000    1.0000    2.0000    2.0000    2.0000
        s[2]    2.0000    2.0000    2.0000    2.0000    2.0000
        s[3]    2.0000    2.0000    2.0000    2.0000    2.0000
In [ ]: