$p(x|w)$ はパラメーター $w$ を持つ $x$ に関する確率密度函数であるとし, $\varphi(w)$ はパラメーター $w$ に関する事前分布(事前確率密度函数)であるとする. サンプル $x=(x_1,\ldots,x_n)$ に対する逆温度 $\beta\geqq 0$ の分配函数は次のように定義される:

$Z^\beta(x) = \int \left(\prod_{i=1}^n p(x_i|w)\right)^\beta\varphi(w)\,dw.$

特に $Z^0(x)=1$ であり, $Z^1(x)$ は通常の分配函数である. このノートでは $\beta\to\infty$ で最尤法が再現されることも数値的に確認する.

逆温度 $\beta$ での事後分布は次のように定義される:

$\Phi^\beta(w|x) = \frac{1}{Z^\beta(x)} \left(\prod_{i=1}^n p(x_i|w)\right)^\beta\varphi(w).$

特に $\Phi^0(w|x)=\varphi(w)$ であり(逆温度 $0$ ＝温度 $\infty$ での事後分布は事前分布と一致), 逆温度 $1$ における事後分布 $\Phi^1(w|x)$ はサンプル $x$ に対する通常の事後分布である. 最尤法が使える状況では $\beta\to\infty$ で事後分布は最尤法の解 $\hat w$ に台を持つデルタ分布に収束する. このノートではそのことを数値的に確認する.

$\beta$ を逆温度と呼ぶ理由は

$H(w|x) = -\sum_{i=1}^n \log p(x_i|w)$

とおいて, 分配函数を

$Z^\beta(x) = \int e^{-\beta H(w|x)}\varphi(w)\,dw.$

と書き直した方が分かり易いかもしれない. これはHamiltonianを $H(w|x)$ の事前分布 $\varphi(w)$ に関するカノニカル分布の分配函数の形になっている. 逆温度 $\beta\to\infty$ の極限すなわち絶対温度→0の極限では, 系の状態は $H(w|x)$ の基底状態 $\hat w$ に凍り付く. ここで $\hat w$ はハミルトニアン $H(w|x)$ の最小値を与える $w$ の値である. $H(w|x)$ は対数尤度函数の $-1$ 倍なので $\hat w$ は最尤法の解(尤度函数の最大値を与えるパラメーター $w$ の値)に一致する.

パッケージの読み込みなど¶

In [ ]:

import PyPlot; plt = PyPlot
using Distributions
using QuadGK
using Mamba

using KernelDensity
function makefunc_pdfkde(X)
    local ik = InterpKDE(kde(X))
    local pdfkde(x) = pdf(ik, x)
    return pdfkde
end
function makefunc_pdfkde(X,Y)
    local ik = InterpKDE(kde((X,Y)))
    local pdfkde(x, y) = pdf(ik, x, y)
    return pdfkde
end

macro sum(f_k, k, itr)
    quote
        begin
            local s = zero(($(esc(k))->$(esc(f_k)))($(esc(itr))[1]))
            for $(esc(k)) in $(esc(itr))
                s += $(esc(f_k))
            end
            s
        end
    end
end

macro prod(f_k, k, itr)
    quote
        begin
            local p = one(($(esc(k))->$(esc(f_k)))($(esc(itr))[1]))
            for $(esc(k)) in $(esc(itr))
                p *= $(esc(f_k))
            end
            p
        end
    end
end

逆温度βでの確率分布の定義¶

Julia言語のMCMCパッケージのMambaで使用する確率モデルの記述はDistributionsパッケージの確率分布によって記述される.

次のセルでは, Distributionsパッケージで定義されている確率分布に対して, 対数確率密度函数(対数尤度函数)の値をβ倍にして返す擬似確率分布を定義する. 逆温度βに依存する擬似確率分布を使って, Mambaで使用する確率モデルを記述すれば, 逆温度βでのベイズ推定が可能になる.

http://mambajl.readthedocs.io/en/latest/mcmc/distributions.html#user-defined-univariate-distributions

In [ ]:

#### Thermal*(d, β) = the distribution d at the inverse temperature β

@everywhere begin
    using Distributions
    import Distributions: logpdf, maximum, minimum
    import Distributions: _logpdf, length, insupport
    import Distributions: pdf, rand, _rand!
    
    mutable struct ThermalCU{T<:Real, D<:ContinuousUnivariateDistribution} <: ContinuousUnivariateDistribution
        d::D
        β::T
    end
    logpdf(d::ThermalCU, x::Real) = d.β*logpdf(d.d, x)
    maximum(d::ThermalCU) = maximum(d.d)
    minimum(d::ThermalCU) = minimum(d.d)

    mutable struct ThermalDU{T<:Real, D<:DiscreteUnivariateDistribution} <: DiscreteUnivariateDistribution
        d::D
        β::T
    end
    logpdf(d::ThermalDU, x::Int) = d.β*logpdf(d.d, x)
    maximum(d::ThermalDU) = maximum(d.d)
    minimum(d::ThermalDU) = minimum(d.d)

    mutable struct ThermalCM{T<:Real, D<:ContinuousMultivariateDistribution} <: ContinuousMultivariateDistribution
        d::D
        β::T
    end
    _logpdf(d::ThermalCM, x::AbstractArray{T,1}) where T = d.β*_logpdf(d.d, x)
    length(d::ThermalCM) = length(d.d)
    insupport(d::ThermalCM, x::AbstractArray{T,1}) where T = insupport(d.d, x)

    mutable struct ThermalDM{T<:Real, D<:DiscreteMultivariateDistribution} <: DiscreteMultivariateDistribution
        d::D
        β::T
    end
    _logpdf(d::ThermalDM, x::AbstractArray{T,1}) where T = d.β*_logpdf(d.d, x)
    length(d::ThermalDM) = length(d.d)
    insupport(d::ThermalDM, x::AbstractArray{T,1}) where T = insupport(d.d, x)
end

### Examples

beta = 0.2

println("-"^70)
@show methods(ThermalCU)
println()
@show dc = ThermalCU(Normal(), beta)
@show logpdf(dc.d, 1.0), logpdf(dc, 1.0)/beta

println("-"^70)
@show methods(ThermalDU)
println()
@show dd = ThermalDU(Poisson(), beta)
@show logpdf(dd.d, 5), logpdf(dd, 5)/beta

println("-"^70)
@show methods(ThermalCM)
println()
@show dcm = ThermalCM(MvNormal([0.0, 0.0], [1.0, 1.0]), beta)
@show _logpdf(dcm.d, [1.0, 2.0]), _logpdf(dcm, [1.0, 2.0])/beta

println("-"^70)
@show methods(ThermalDM)
println()
@show ddm = ThermalDM(Multinomial(10, [0.1, 0.3, 0.6]), beta)
@show _logpdf(ddm.d, [1,3,6]), _logpdf(ddm, [1,3,6])/beta
;

MCMCの結果のまとめの表示函数¶

In [ ]:

## Summary
function showsummary(sim;
        sortkeys=true, figsize_t=(8, 3), figsize_c=(8, 3.5), maxfigsize=8.0,
        show_describe=false, show_gelman=false, plot_trace=true, plot_contour=false, plot_chain=true)
    ## Summary of MCMC
    if show_describe
        println("\n========== Summary:\n")
        display(describe(sim))
    end

    # Convergence Diagnostics
    if show_gelman && length(sim.chains) ≥ 2 
       println("========== Gelman Diagnostic:")
       show(gelmandiag(sim))
    end

    ## Plot
    sleep(0.1)
    if plot_trace
        #draw(plot(sim), fmt=:png, width=10inch, height=3.5inch, nrow=1, ncol=2, ask=false)
        pyplot_trace(sim, sortkeys=sortkeys, figsize=figsize_t)
    end
    if plot_contour
        #draw(plot(sim, :contour), fmt=:png, width=10inch, height=4.5inch, nrow=1, ncol=2, ask=false)
        pyplot_contour(sim, sortkeys=sortkeys, figsize=figsize_c)
    end
    if plot_chain
        plotchain(sim, maxfigsize=maxfigsize)
    end
end

## plot traces
function pyplot_trace(sim; sortkeys = false, figsize = (8, 3))
    if sortkeys
        keys_sim = sort(keys(sim))
    else
        keys_sim = keys(sim)
    end
    for var in keys_sim
        plt.figure(figsize=figsize)
        
        plt.subplot(1,2,1)
        for k in sim.chains
            plt.plot(sim.range, sim[:,var,:].value[:,1,k], label="$k", lw=0.2, alpha=0.8)
        end
        plt.xlabel("iterations")
        plt.grid(ls=":")
        #plt.legend(loc="upper right")
        plt.title("trace of $var", fontsize=10)
        
        plt.subplot(1,2,2)
        local xmin = quantile(vec(sim[:,var,:].value), 0.005)
        local xmax = quantile(vec(sim[:,var,:].value), 0.995)
        for k in sim.chains
            local chain = sim[:,var,:].value[:,1,k]
            local pdfkde = makefunc_pdfkde(chain)
            local x = linspace(xmin, xmax, 201)
            plt.plot(x, pdfkde.(x), label="$k", lw=0.8, alpha=0.8)
        end
        plt.xlabel("$var")
        plt.grid(ls=":")
        plt.title("empirical posterior pdf of $var", fontsize=10)

        plt.tight_layout()
    end
end

# plot contours
function pyplot_contour(sim; sortkeys = true, figsize = (8, 3.5))
    if sortkeys
        keys_sim = sort(keys(sim))
    else
        keys_sim = keys(sim)
    end
    local m = 0
    local K = length(keys_sim)
    for i in 1:K
        for j in i+1:K
            m += 1
            mod1(m,2) == 1 && plt.figure(figsize=figsize)
            plt.subplot(1,2, mod1(m,2))

            local u = keys_sim[i]
            local v = keys_sim[j]
            local X = vec(sim[:,u,:].value)
            local Y = vec(sim[:,v,:].value)
            local pdfkde = makefunc_pdfkde(X,Y)
            local xmin = quantile(X, 0.005)
            local xmax = quantile(X, 0.995)
            local ymin = quantile(Y, 0.005)
            local ymax = quantile(Y, 0.995)
            local x = linspace(xmin, xmax, 200)
            local y = linspace(ymin, ymax, 200)
            
            plt.pcolormesh(x', y, pdfkde.(x',y), cmap="CMRmap")
            plt.colorbar()
            plt.grid(ls=":")
            plt.xlabel(u)
            plt.ylabel(v)
            plt.title("posterior of ($u, $v)", fontsize=10)

            mod1(m,2) == 2 && plt.tight_layout()

            if 2*m == K*(K-1) && mod1(m,2) == 1
                plt.subplot(1,2,2)
                
                plt.pcolormesh(y', x, pdfkde.(x,y'), cmap="CMRmap")
                plt.colorbar()
                plt.grid(ls=":")
                plt.xlabel(v)
                plt.ylabel(u)
                plt.title("posterior of ($v, $u)", fontsize=10)

                plt.tight_layout()
            end
        end
    end
end

# plot chains
function plotchain(chain, paramname; maxfigsize=8.0)
    local K = size(chain,2)
    local m = 0
    local fs = min(maxfigsize, 3.0*K)
    plt.figure(figsize=(1.1*fs,fs))
    for j in 1:K
        for i in 1:K
            m += 1
            fig = plt.subplot(K,K,m)
            fig[:set_facecolor]("k")
            if i == j
                x0 = minimum(@view chain[:,i])
                x1 = maximum(@view chain[:,i])
                x = linspace(x0, x1, 200)
                ik = InterpKDE(kde(chain[:,i]))
                plt.plot(x, pdf.(ik,x), color="pink", lw=1)
                xmin = quantile((@view chain[:,i]), 0.005)
                xmax = quantile((@view chain[:,i]), 0.995)
                plt.xlim(xmin, xmax)
                plt.grid(ls=":")
                plt.xticks(fontsize=6)
                plt.yticks(fontsize=6)
                j == K && plt.xlabel("$(paramname[i])", fontsize=10)
                i == 1 && plt.ylabel("$(paramname[i])", fontsize=10)
            else
                ik = InterpKDE(kde((chain[:,i], chain[:,j])))
                plt.scatter(chain[:,i], chain[:,j], c=pdf.(ik, chain[:,i], chain[:,j]), s=2, cmap="CMRmap")
                plt.grid(ls=":")
                xmin = quantile((@view chain[:,i]), 0.005)
                xmax = quantile((@view chain[:,i]), 0.995)
                ymin = quantile((@view chain[:,j]), 0.005)
                ymax = quantile((@view chain[:,j]), 0.995)
                plt.xlim(xmin, xmax)
                plt.ylim(ymin, ymax)
                plt.xticks(fontsize=6)
                plt.yticks(fontsize=6)
                j == K && plt.xlabel("$(paramname[i])", fontsize=10)
                i == 1 && plt.ylabel("$(paramname[j])", fontsize=10)
            end
        end
    end
    plt.tight_layout()
end

sim2chain(sim, paramname) = hcat((vec(sim[:,name,:].value) for name in paramname)...)

function plotchain(sim::ModelChains, paramname; maxfigsize=8.0)
    plotchain(sim2chain(sim, paramname), paramname, maxfigsize=maxfigsize)
end

function plotchain(sim::ModelChains; maxfigsize=8.0)
    local paramname = sort(keys(sim))
    plotchain(sim2chain(sim, paramname), paramname, maxfigsize=maxfigsize)
end

分散を1に固定した正規分布モデル normal1¶

In [ ]:

@everywhere normal1(mu::Real) = Normal(mu, 1.0)
@everywhere thermal_normal1(mu, beta) = ThermalCU(normal1(mu), beta)

normal1(w::AbstractVector) = normal1(w[1])
unlink_normal1(w) = w
link_normal1(z)   = z

pdf_normal1(mu, x) = pdf(normal1(w), x)
function rand_prior_normal(M)
    local mu = randn(M)
    local sigma = abs.(randn(M))
    return [mu';sigma']
end

function sample2model_normal1(Y;
        beta = 1.0,
        thermal_dist_model = thermal_normal1,
        mu0    = 0.0,
        prior_mu = Normal(0.0, 1.0),
        chains = 2,
    )
    local data = Dict{Symbol, Any}(
        :Y => Y,
        :n => length(Y),
        :mu0 => mu0,
        :beta => beta
    )
    
    local model = Model(
        y = Stochastic(1, mu -> thermal_dist_model(mu, data[:beta]), false),
        mu = Stochastic(() -> prior_mu, true),
    )
    local scheme = [
        NUTS([:mu])
    ]
    setsamplers!(model, scheme)
    
    local inits = [
        Dict{Symbol, Any}(
            :y => data[:Y],
            :mu => data[:mu0],
        )
        for k in 1:chains
    ]
    return model, data, inits
end

srand(12345)
dist_true = Normal(0.5,1.0)
n = 2^7
Y = rand(dist_true, n)
beta = 1.0

println("========== Sample")
@show dist_true
@show n
@show mean(Y)
@show std(Y)
println("========== beta")
@show beta
println("="^20)

iterations = 4200
burnin = 200
thin = 1
chains = 4
model, data, inits = sample2model_normal1(Y, beta=beta, chains=chains)
@time sim_normal1_1 = mcmc(model, data, inits, iterations, burnin=burnin, thin=thin, chains=chains, verbose=false)
sim = sim_normal1_1
showsummary(sim, show_describe=true, show_gelman=false, plot_trace=true, plot_chain=true)

In [ ]:

@everywhere normal1(mu::Real) = Normal(mu, 1.0)
@everywhere thermal_normal1(mu, beta) = ThermalCU(normal1(mu), beta)

normal1(w::AbstractVector) = normal1(w[1])
unlink_normal1(w) = w
link_normal1(z)   = z

pdf_normal1(mu, x) = pdf(normal1(w), x)
function rand_prior_normal(M)
    local mu = randn(M)
    local sigma = abs.(randn(M))
    return [mu';sigma']
end

function sample2model_normal1(Y;
        beta = 1.0,
        thermal_dist_model = thermal_normal1,
        mu0    = 0.0,
        prior_mu = Normal(0.0, 1.0),
        chains = 2,
    )
    local data = Dict{Symbol, Any}(
        :Y => Y,
        :n => length(Y),
        :mu0 => mu0,
        :beta => beta
    )
    
    local model = Model(
        y = Stochastic(1, mu -> thermal_dist_model(mu, data[:beta]), false),
        mu = Stochastic(() -> prior_mu, true),
    )
    local scheme = [
        NUTS([:mu])
    ]
    setsamplers!(model, scheme)
    
    local inits = [
        Dict{Symbol, Any}(
            :y => data[:Y],
            :mu => data[:mu0],
        )
        for k in 1:chains
    ]
    return model, data, inits
end

srand(12345)
dist_true = Normal(0.5,1.0)
n = 2^7
Y = rand(dist_true, n)
beta = 0.0

println("========== Sample")
@show dist_true
@show n
@show mean(Y)
@show std(Y)
println("========== beta")
@show beta
println("="^20)

iterations = 4200
burnin = 200
thin = 1
chains = 4
model, data, inits = sample2model_normal1(Y, beta=beta, chains=chains)
@time sim_normal1_0 = mcmc(model, data, inits, iterations, burnin=burnin, thin=thin, chains=chains, verbose=false)
sim = sim_normal1_0
showsummary(sim, show_describe=true, show_gelman=false, plot_trace=true, plot_chain=true)

In [ ]:

@everywhere normal1(mu::Real) = Normal(mu, 1.0)
@everywhere thermal_normal1(mu, beta) = ThermalCU(normal1(mu), beta)

normal1(w::AbstractVector) = normal1(w[1])
unlink_normal1(w) = w
link_normal1(z)   = z

pdf_normal1(mu, x) = pdf(normal1(w), x)
function rand_prior_normal(M)
    local mu = randn(M)
    local sigma = abs.(randn(M))
    return [mu';sigma']
end

function sample2model_normal1(Y;
        beta = 1.0,
        thermal_dist_model = thermal_normal1,
        mu0    = 0.0,
        prior_mu = Normal(0.0, 1.0),
        chains = 2,
    )
    local data = Dict{Symbol, Any}(
        :Y => Y,
        :n => length(Y),
        :mu0 => mu0,
        :beta => beta
    )
    
    local model = Model(
        y = Stochastic(1, mu -> thermal_dist_model(mu, data[:beta]), false),
        mu = Stochastic(() -> prior_mu, true),
    )
    local scheme = [
        NUTS([:mu])
    ]
    setsamplers!(model, scheme)
    
    local inits = [
        Dict{Symbol, Any}(
            :y => data[:Y],
            :mu => data[:mu0],
        )
        for k in 1:chains
    ]
    return model, data, inits
end

srand(12345)
dist_true = Normal(0.5,1.0)
n = 2^7
Y = rand(dist_true, n)
beta = 10000.0

println("========== Sample")
@show dist_true
@show n
@show mean(Y)
@show std(Y)
println("========== beta")
@show beta
println("="^20)

iterations = 4200
burnin = 200
thin = 1
chains = 4
model, data, inits = sample2model_normal1(Y, beta=beta, chains=chains)
@time sim_normal1_10000 = mcmc(model, data, inits, iterations, burnin=burnin, thin=thin, chains=chains, verbose=false)
sim = sim_normal1_10000
showsummary(sim, show_describe=true, show_gelman=false, plot_trace=true, plot_chain=true)

逆温度を $\beta=10000.0$ と大きな値にすると, 事後分布はサンプルの平均 $\mu=0.508$ に台を持つデルタ分布に近い分布になる. 一般に逆温度 $\beta\to\infty$ の極限で事後分布は最尤法の解に台を持つデルタ分布に収束する.

normal1モデルと同等のモデルを多変量正規分布を使って表現した場合¶

In [ ]:

function sample2model_mvnormal1(Y;
        beta = 1.0,
        thermal_dist_model = Flat,
        mu0    = 0.0,
        prior_mu = Normal(0.0, 1.0),
        chains = 2,
    )
    local data = Dict{Symbol, Any}(
        :Y => Y,
        :n => length(Y),
        :mu0 => mu0,
        :beta => beta
    )
    
    local model = Model(
        y = Stochastic(1, mu -> ThermalCM(MvNormal(fill(mu, data[:n]), ones(data[:n])), data[:beta]), false),
        mu = Stochastic(() -> prior_mu, true),
    )
    local scheme = [
        NUTS([:mu])
    ]
    setsamplers!(model, scheme)
    
    local inits = [
        Dict{Symbol, Any}(
            :y => data[:Y],
            :mu => data[:mu0],
        )
        for k in 1:chains
    ]
    return model, data, inits
end

srand(12345)
dist_true = Normal(0.5,1.0)
n = 2^7
Y = rand(dist_true, n)
beta = 10000.0

println("========== Sample")
@show dist_true
@show n
@show mean(Y)
@show std(Y)
println("========== beta")
@show beta
println("="^20)

iterations = 4200
burnin = 200
thin = 1
chains = 4
model, data, inits = sample2model_mvnormal1(Y, beta=beta, chains=chains)
@time sim_mvnormal_10000 = mcmc(model, data, inits, iterations, burnin=burnin, thin=thin, chains=chains, verbose=false)
sim = sim_mvnormal_10000
showsummary(sim, show_describe=true, show_gelman=false, plot_trace=true, plot_chain=true)

分散1の2つの正規分布の混合モデル mixnormal¶

In [ ]:

@everywhere mixnormal(a,b,c) = MixtureModel(Normal[Normal(b, 1.0), Normal(c, 1.0)], [1.0-a, a])
@everywhere thermal_mixnormal(a, b, c, beta) = ThermalCU(mixnormal(a,b,c), beta)

mixnormal(w::AbstractVector) = mixnormal(w[1], w[2], w[3])
unlink_mixnormal(w) = [logit(w[1]), w[2], w[3]]     # unlink_mixnormal : (0,1)×R^2 -> R^3
link_mixnormal(z)   = [invlogit(z[1]), z[2], z[3]]  #   link_mixnormal : R^3 → (0,1)×R^2

pdf_mixnormal(w,x) = (1-w[1])*pdf(normal1(w[2]),x) + w[1]*pdf(normal1(w[3]),x)
function rand_prior_mixnormal(M)
    local a = rand(M)
    local b = randn(M)
    local c = randn(M)
    return [a';b';c']
end

function sample2model_mixnormal(Y;
        beta = 1.0,
        dist_model = thermal_mixnormal,
        a0 = 0.5,
        b0 = 0.0,
        c0 = 0.0,
        prior_a = Uniform(0.0, 1.0),
        prior_b = Normal(0.0, 1.0),
        prior_c = Normal(0.0, 1.0),
        chains = 2
    )
    local data = Dict{Symbol, Any}(
        :Y => Y,
        :n => length(Y),
        :a0 => a0,
        :b0 => b0,
        :c0 => c0,
        :beta => beta,
    )
    
    local model = Model(
        y = Stochastic(1, (a, b, c) -> dist_model(a, b, c, data[:beta]), false),
        a = Stochastic(() -> prior_a, true),
        b = Stochastic(() -> prior_b, true),
        c = Stochastic(() -> prior_b, true),
    )
    local scheme = [
        NUTS([:a, :b, :c])
    ]
    setsamplers!(model, scheme)
    
    local inits = [
        Dict{Symbol, Any}(
            :y => data[:Y],
            :a => data[:a0],
            :b => data[:b0],
            :c => data[:c0],
        )
        for k in 1:chains
    ]
    return model, data, inits
end

srand(12345)
dist_true = Normal(0.0,1.0)
n = 2^7
Y = rand(dist_true, n)
beta = 1.0

println("========== Sample")
@show dist_true
@show n
@show mean(Y)
@show std(Y)
println("========== beta")
@show beta
println("="^20)

iterations = 1200
burnin = 200
thin = 1
chains = 4
model, data, inits = sample2model_mixnormal(Y, beta=beta, chains=chains)
@time sim_mixnormal_1 = mcmc(model, data, inits, iterations, burnin=burnin, thin=thin, chains=chains, verbose=false)
sim = sim_mixnormal_1
showsummary(sim, show_describe=true, show_gelman=false, plot_trace=true, plot_chain=true)

擬似確率分布を直接定義することによって高速化¶

上の推定では mixnormal モデルを Distributions パッケージの MixtureModel を使って定義した.

しかし, 以下のように mixnormal モデルを直接定義すると計算速度が10倍程度速くなる!

In [ ]:

@everywhere mutable struct ThermalMixNormal1{T<:Real} <:ContinuousUnivariateDistribution
    a::T
    b::T
    c::T
    β::T
end
@everywhere ThermalMixNormal1(a::Real, b::Real, c::Real, β::Real) =
    ThermalMixNormal1(float(a), float(b), float(c), float(β))
@everywhere function Distributions.pdf(d::ThermalMixNormal1, x::Real)
    local p = (1.0 - d.a)*exp(-0.5*(x-d.b)^2) + d.a*exp(-0.5*(x-d.c)^2) 
    p /= sqrt(2*pi)
    return p^d.β
end
@everywhere function Distributions.logpdf(d::ThermalMixNormal1, x::Real)
    local l = log((1.0 - d.a)*exp(-0.5*(x-d.b)^2) + d.a*exp(-0.5*(x-d.c)^2))
    l -= 0.5*log(2*pi)
    return d.β*l
end
@everywhere Distributions.maximum(d::ThermalMixNormal1) = Inf
@everywhere Distributions.minimum(d::ThermalMixNormal1) = -Inf

mixnormal1(w::AbstractVector) = ThermalMixNormal1(w[1], w[2], w[3], 1.0)
unlink_mixnormal1(w) = [logit(w[1]), w[2], w[3]]     # unlink_mixnormal : (0,1)×R^2 -> R^3
link_mixnormal1(z)   = [invlogit(z[1]), z[2], z[3]]  #   link_mixnormal : R^3 → (0,1)×R^2

pdf_mixnormal(w,x) = pdf(ThermalMixNormal1(w[1], w[2], w[3], 1.0), x)
function rand_prior_mixnormal(M)
    local a = rand(M)
    local b = randn(M)
    local c = randn(M)
    return [a';b';c']
end

function sample2model_mixnormal1(Y;
        beta = 1.0,
        dist_model = ThermalMixNormal1,
        a0 = 0.5,
        b0 = 0.0,
        c0 = 0.0,
        prior_a = Uniform(0.0, 1.0),
        prior_b = Normal(0.0, 1.0),
        prior_c = Normal(0.0, 1.0),
        chains = 2
    )
    local data = Dict{Symbol, Any}(
        :Y => Y,
        :n => length(Y),
        :a0 => a0,
        :b0 => b0,
        :c0 => c0,
        :beta => beta,
    )
    
    local model = Model(
        y = Stochastic(1, (a, b, c) -> dist_model(a, b, c, data[:beta]), false),
        a = Stochastic(() -> prior_a, true),
        b = Stochastic(() -> prior_b, true),
        c = Stochastic(() -> prior_b, true),
    )
    local scheme = [
        NUTS([:a, :b, :c])
    ]
    setsamplers!(model, scheme)
    
    local inits = [
        Dict{Symbol, Any}(
            :y => data[:Y],
            :a => data[:a0],
            :b => data[:b0],
            :c => data[:c0],
        )
        for k in 1:chains
    ]
    return model, data, inits
end

srand(12345)
dist_true = Normal(0.0,1.0)
n = 2^7
Y = rand(dist_true, n)
beta = 1.0

println("========== Sample")
@show dist_true
@show n
@show mean(Y)
@show std(Y)
println("========== beta")
@show beta
println("="^20)

iterations = 1200
burnin = 200
thin = 1
chains = 4
model, data, inits = sample2model_mixnormal1(Y, beta=beta, chains=chains)
@time sim_mixnormal_1 = mcmc(model, data, inits, iterations, burnin=burnin, thin=thin, chains=chains, verbose=false)
sim = sim_mixnormal_1
showsummary(sim, show_describe=true, show_gelman=false, plot_trace=true, plot_chain=true)

高速化に成功したので、iterations = 4200 でMCMCを回してみる。¶

In [ ]:

@everywhere mutable struct ThermalMixNormal1{T<:Real} <:ContinuousUnivariateDistribution
    a::T
    b::T
    c::T
    β::T
end
@everywhere ThermalMixNormal1(a::Real, b::Real, c::Real, β::Real) =
    ThermalMixNormal1(float(a), float(b), float(c), float(β))
@everywhere function Distributions.pdf(d::ThermalMixNormal1, x::Real)
    local p = (1.0 - d.a)*exp(-0.5*(x-d.b)^2) + d.a*exp(-0.5*(x-d.c)^2) 
    p /= sqrt(2*pi)
    return p^d.β
end
@everywhere function Distributions.logpdf(d::ThermalMixNormal1, x::Real)
    local l = log((1.0 - d.a)*exp(-0.5*(x-d.b)^2) + d.a*exp(-0.5*(x-d.c)^2))
    l -= 0.5*log(2*pi)
    return d.β*l
end
@everywhere Distributions.maximum(d::ThermalMixNormal1) = Inf
@everywhere Distributions.minimum(d::ThermalMixNormal1) = -Inf

mixnormal1(w::AbstractVector) = ThermalMixNormal1(w[1], w[2], w[3], 1.0)
unlink_mixnormal1(w) = [logit(w[1]), w[2], w[3]]     # unlink_mixnormal : (0,1)×R^2 -> R^3
link_mixnormal1(z)   = [invlogit(z[1]), z[2], z[3]]  #   link_mixnormal : R^3 → (0,1)×R^2

pdf_mixnormal(w,x) = pdf(ThermalMixNormal1(w[1], w[2], w[3], 1.0), x)
function rand_prior_mixnormal(M)
    local a = rand(M)
    local b = randn(M)
    local c = randn(M)
    return [a';b';c']
end

function sample2model_mixnormal1(Y;
        beta = 1.0,
        dist_model = ThermalMixNormal1,
        a0 = 0.5,
        b0 = 0.0,
        c0 = 0.0,
        prior_a = Uniform(0.0, 1.0),
        prior_b = Normal(0.0, 1.0),
        prior_c = Normal(0.0, 1.0),
        chains = 2
    )
    local data = Dict{Symbol, Any}(
        :Y => Y,
        :n => length(Y),
        :a0 => a0,
        :b0 => b0,
        :c0 => c0,
        :beta => beta,
    )
    
    local model = Model(
        y = Stochastic(1, (a, b, c) -> dist_model(a, b, c, data[:beta]), false),
        a = Stochastic(() -> prior_a, true),
        b = Stochastic(() -> prior_b, true),
        c = Stochastic(() -> prior_b, true),
    )
    local scheme = [
        NUTS([:a, :b, :c])
    ]
    setsamplers!(model, scheme)
    
    local inits = [
        Dict{Symbol, Any}(
            :y => data[:Y],
            :a => data[:a0],
            :b => data[:b0],
            :c => data[:c0],
        )
        for k in 1:chains
    ]
    return model, data, inits
end

srand(12345)
dist_true = Normal(0.0,1.0)
n = 2^7
Y = rand(dist_true, n)
beta = 1.0

println("========== Sample")
@show dist_true
@show n
@show mean(Y)
@show std(Y)
println("========== beta")
@show beta
println("="^20)

iterations = 4200
burnin = 200
thin = 1
chains = 4
model, data, inits = sample2model_mixnormal1(Y, beta=beta, chains=chains)
@time sim_mixnormal_1 = mcmc(model, data, inits, iterations, burnin=burnin, thin=thin, chains=chains, verbose=false)
sim = sim_mixnormal_1
showsummary(sim, show_describe=true, show_gelman=false, plot_trace=true, plot_chain=true)

In [ ]:

@everywhere mutable struct ThermalMixNormal1{T<:Real} <:ContinuousUnivariateDistribution
    a::T
    b::T
    c::T
    β::T
end
@everywhere ThermalMixNormal1(a::Real, b::Real, c::Real, β::Real) =
    ThermalMixNormal1(float(a), float(b), float(c), float(β))
@everywhere function Distributions.pdf(d::ThermalMixNormal1, x::Real)
    local p = (1.0 - d.a)*exp(-0.5*(x-d.b)^2) + d.a*exp(-0.5*(x-d.c)^2) 
    p /= sqrt(2*pi)
    return p^d.β
end
@everywhere function Distributions.logpdf(d::ThermalMixNormal1, x::Real)
    local l = log((1.0 - d.a)*exp(-0.5*(x-d.b)^2) + d.a*exp(-0.5*(x-d.c)^2))
    l -= 0.5*log(2*pi)
    return d.β*l
end
@everywhere Distributions.maximum(d::ThermalMixNormal1) = Inf
@everywhere Distributions.minimum(d::ThermalMixNormal1) = -Inf

mixnormal1(w::AbstractVector) = ThermalMixNormal1(w[1], w[2], w[3], 1.0)
unlink_mixnormal1(w) = [logit(w[1]), w[2], w[3]]     # unlink_mixnormal : (0,1)×R^2 -> R^3
link_mixnormal1(z)   = [invlogit(z[1]), z[2], z[3]]  #   link_mixnormal : R^3 → (0,1)×R^2

pdf_mixnormal(w,x) = pdf(ThermalMixNormal1(w[1], w[2], w[3], 1.0), x)
function rand_prior_mixnormal(M)
    local a = rand(M)
    local b = randn(M)
    local c = randn(M)
    return [a';b';c']
end

function sample2model_mixnormal1(Y;
        beta = 1.0,
        dist_model = ThermalMixNormal1,
        a0 = 0.5,
        b0 = 0.0,
        c0 = 0.0,
        prior_a = Uniform(0.0, 1.0),
        prior_b = Normal(0.0, 1.0),
        prior_c = Normal(0.0, 1.0),
        chains = 2
    )
    local data = Dict{Symbol, Any}(
        :Y => Y,
        :n => length(Y),
        :a0 => a0,
        :b0 => b0,
        :c0 => c0,
        :beta => beta,
    )
    
    local model = Model(
        y = Stochastic(1, (a, b, c) -> dist_model(a, b, c, data[:beta]), false),
        a = Stochastic(() -> prior_a, true),
        b = Stochastic(() -> prior_b, true),
        c = Stochastic(() -> prior_b, true),
    )
    local scheme = [
        NUTS([:a, :b, :c])
    ]
    setsamplers!(model, scheme)
    
    local inits = [
        Dict{Symbol, Any}(
            :y => data[:Y],
            :a => data[:a0],
            :b => data[:b0],
            :c => data[:c0],
        )
        for k in 1:chains
    ]
    return model, data, inits
end

srand(12345)
dist_true = Normal(0.0,1.0)
n = 2^7
Y = rand(dist_true, n)
beta = 1.0/log(n)

println("========== Sample")
@show dist_true
@show n
@show mean(Y)
@show std(Y)
println("========== beta")
@show beta
println("="^20)

iterations = 4200
burnin = 200
thin = 1
chains = 4
model, data, inits = sample2model_mixnormal1(Y, beta=beta, chains=chains)
@time sim_mixnormal1_invlogn = mcmc(model, data, inits, iterations, burnin=burnin, thin=thin, chains=chains, verbose=false)
sim = sim_mixnormal1_invlogn
showsummary(sim, show_describe=true, show_gelman=false, plot_trace=true, plot_chain=true)

In [ ]:

@everywhere mutable struct ThermalMixNormal1{T<:Real} <:ContinuousUnivariateDistribution
    a::T
    b::T
    c::T
    β::T
end
@everywhere ThermalMixNormal1(a::Real, b::Real, c::Real, β::Real) =
    ThermalMixNormal1(float(a), float(b), float(c), float(β))
@everywhere function Distributions.pdf(d::ThermalMixNormal1, x::Real)
    local p = (1.0 - d.a)*exp(-0.5*(x-d.b)^2) + d.a*exp(-0.5*(x-d.c)^2) 
    p /= sqrt(2*pi)
    return p^d.β
end
@everywhere function Distributions.logpdf(d::ThermalMixNormal1, x::Real)
    local l = log((1.0 - d.a)*exp(-0.5*(x-d.b)^2) + d.a*exp(-0.5*(x-d.c)^2))
    l -= 0.5*log(2*pi)
    return d.β*l
end
@everywhere Distributions.maximum(d::ThermalMixNormal1) = Inf
@everywhere Distributions.minimum(d::ThermalMixNormal1) = -Inf

mixnormal1(w::AbstractVector) = ThermalMixNormal1(w[1], w[2], w[3], 1.0)
unlink_mixnormal1(w) = [logit(w[1]), w[2], w[3]]     # unlink_mixnormal : (0,1)×R^2 -> R^3
link_mixnormal1(z)   = [invlogit(z[1]), z[2], z[3]]  #   link_mixnormal : R^3 → (0,1)×R^2

pdf_mixnormal(w,x) = pdf(ThermalMixNormal1(w[1], w[2], w[3], 1.0), x)
function rand_prior_mixnormal(M)
    local a = rand(M)
    local b = randn(M)
    local c = randn(M)
    return [a';b';c']
end

function sample2model_mixnormal1(Y;
        beta = 1.0,
        dist_model = ThermalMixNormal1,
        a0 = 0.5,
        b0 = 0.0,
        c0 = 0.0,
        prior_a = Uniform(0.0, 1.0),
        prior_b = Normal(0.0, 1.0),
        prior_c = Normal(0.0, 1.0),
        chains = 2
    )
    local data = Dict{Symbol, Any}(
        :Y => Y,
        :n => length(Y),
        :a0 => a0,
        :b0 => b0,
        :c0 => c0,
        :beta => beta,
    )
    
    local model = Model(
        y = Stochastic(1, (a, b, c) -> dist_model(a, b, c, data[:beta]), false),
        a = Stochastic(() -> prior_a, true),
        b = Stochastic(() -> prior_b, true),
        c = Stochastic(() -> prior_b, true),
    )
    local scheme = [
        NUTS([:a, :b, :c])
    ]
    setsamplers!(model, scheme)
    
    local inits = [
        Dict{Symbol, Any}(
            :y => data[:Y],
            :a => data[:a0],
            :b => data[:b0],
            :c => data[:c0],
        )
        for k in 1:chains
    ]
    return model, data, inits
end

srand(12345)
dist_true = Normal(0.0,1.0)
n = 2^7
Y = rand(dist_true, n)
beta = 100.0

println("========== Sample")
@show dist_true
@show n
@show mean(Y)
@show std(Y)
println("========== beta")
@show beta
println("="^20)

iterations = 4200
burnin = 200
thin = 1
chains = 4
model, data, inits = sample2model_mixnormal1(Y, beta=beta, chains=chains)
@time sim_mixnormal1_100 = mcmc(model, data, inits, iterations, burnin=burnin, thin=thin, chains=chains, verbose=false)
sim = sim_mixnormal1_100
showsummary(sim, show_describe=true, show_gelman=false, plot_trace=true, plot_chain=true)

サンプルを二つ山の分布で生成してみる¶

In [ ]:

@everywhere mutable struct ThermalMixNormal1{T<:Real} <:ContinuousUnivariateDistribution
    a::T
    b::T
    c::T
    β::T
end
@everywhere ThermalMixNormal1(a::Real, b::Real, c::Real, β::Real) =
    ThermalMixNormal1(float(a), float(b), float(c), float(β))
@everywhere function Distributions.pdf(d::ThermalMixNormal1, x::Real)
    local p = (1.0 - d.a)*exp(-0.5*(x-d.b)^2) + d.a*exp(-0.5*(x-d.c)^2) 
    p /= sqrt(2*pi)
    return p^d.β
end
@everywhere function Distributions.logpdf(d::ThermalMixNormal1, x::Real)
    local l = log((1.0 - d.a)*exp(-0.5*(x-d.b)^2) + d.a*exp(-0.5*(x-d.c)^2))
    l -= 0.5*log(2*pi)
    return d.β*l
end
@everywhere Distributions.maximum(d::ThermalMixNormal1) = Inf
@everywhere Distributions.minimum(d::ThermalMixNormal1) = -Inf

mixnormal1(w::AbstractVector) = ThermalMixNormal1(w[1], w[2], w[3], 1.0)
unlink_mixnormal1(w) = [logit(w[1]), w[2], w[3]]     # unlink_mixnormal : (0,1)×R^2 -> R^3
link_mixnormal1(z)   = [invlogit(z[1]), z[2], z[3]]  #   link_mixnormal : R^3 → (0,1)×R^2

pdf_mixnormal(w,x) = pdf(ThermalMixNormal1(w[1], w[2], w[3], 1.0), x)
function rand_prior_mixnormal(M)
    local a = rand(M)
    local b = randn(M)
    local c = randn(M)
    return [a';b';c']
end

function sample2model_mixnormal1(Y;
        beta = 1.0,
        dist_model = ThermalMixNormal1,
        a0 = 0.5,
        b0 = 0.0,
        c0 = 0.0,
        prior_a = Uniform(0.0, 1.0),
        prior_b = Normal(0.0, 2.0),
        prior_c = Normal(0.0, 2.0),
        chains = 2
    )
    local data = Dict{Symbol, Any}(
        :Y => Y,
        :n => length(Y),
        :a0 => a0,
        :b0 => b0,
        :c0 => c0,
        :beta => beta,
    )
    
    local model = Model(
        y = Stochastic(1, (a, b, c) -> dist_model(a, b, c, data[:beta]), false),
        a = Stochastic(() -> prior_a, true),
        b = Stochastic(() -> prior_b, true),
        c = Stochastic(() -> prior_b, true),
    )
    local scheme = [
        NUTS([:a, :b, :c])
    ]
    setsamplers!(model, scheme)
    
    local inits = [
        Dict{Symbol, Any}(
            :y => data[:Y],
            :a => data[:a0],
            :b => data[:b0],
            :c => data[:c0],
        )
        for k in 1:chains
    ]
    return model, data, inits
end

srand(12345)
dist_true = mixnormal(0.3, -0.5, 0.5)
n = 2^7
Y = rand(dist_true, n)
beta = 1.0

println("========== Sample")
@show dist_true
@show n
@show mean(Y)
@show std(Y)
println("========== beta")
@show beta
println("="^20)

iterations = 4200
burnin = 200
thin = 1
chains = 8
model, data, inits = sample2model_mixnormal1(Y, beta=beta, chains=chains)
@time sim = mcmc(model, data, inits, iterations, burnin=burnin, thin=thin, chains=chains, verbose=false)
showsummary(sim, show_describe=true, show_gelman=false, plot_trace=true, plot_chain=true)

サンプルを三つ山の分布で生成してみる¶

In [ ]:

@everywhere mutable struct ThermalMixNormal1{T<:Real} <:ContinuousUnivariateDistribution
    a::T
    b::T
    c::T
    β::T
end
@everywhere ThermalMixNormal1(a::Real, b::Real, c::Real, β::Real) =
    ThermalMixNormal1(float(a), float(b), float(c), float(β))
@everywhere function Distributions.pdf(d::ThermalMixNormal1, x::Real)
    local p = (1.0 - d.a)*exp(-0.5*(x-d.b)^2) + d.a*exp(-0.5*(x-d.c)^2) 
    p /= sqrt(2*pi)
    return p^d.β
end
@everywhere function Distributions.logpdf(d::ThermalMixNormal1, x::Real)
    local l = log((1.0 - d.a)*exp(-0.5*(x-d.b)^2) + d.a*exp(-0.5*(x-d.c)^2))
    l -= 0.5*log(2*pi)
    return d.β*l
end
@everywhere Distributions.maximum(d::ThermalMixNormal1) = Inf
@everywhere Distributions.minimum(d::ThermalMixNormal1) = -Inf

mixnormal1(w::AbstractVector) = ThermalMixNormal1(w[1], w[2], w[3], 1.0)
unlink_mixnormal1(w) = [logit(w[1]), w[2], w[3]]     # unlink_mixnormal : (0,1)×R^2 -> R^3
link_mixnormal1(z)   = [invlogit(z[1]), z[2], z[3]]  #   link_mixnormal : R^3 → (0,1)×R^2

pdf_mixnormal(w,x) = pdf(ThermalMixNormal1(w[1], w[2], w[3], 1.0), x)
function rand_prior_mixnormal(M)
    local a = rand(M)
    local b = randn(M)
    local c = randn(M)
    return [a';b';c']
end

function sample2model_mixnormal1(Y;
        beta = 1.0,
        dist_model = ThermalMixNormal1,
        a0 = 0.5,
        b0 = 0.0,
        c0 = 0.0,
        prior_a = Uniform(0.0, 1.0),
        prior_b = Normal(0.0, 2.0),
        prior_c = Normal(0.0, 2.0),
        chains = 2
    )
    local data = Dict{Symbol, Any}(
        :Y => Y,
        :n => length(Y),
        :a0 => a0,
        :b0 => b0,
        :c0 => c0,
        :beta => beta,
    )
    
    local model = Model(
        y = Stochastic(1, (a, b, c) -> dist_model(a, b, c, data[:beta]), false),
        a = Stochastic(() -> prior_a, true),
        b = Stochastic(() -> prior_b, true),
        c = Stochastic(() -> prior_b, true),
    )
    local scheme = [
        NUTS([:a, :b, :c])
    ]
    setsamplers!(model, scheme)
    
    local inits = [
        Dict{Symbol, Any}(
            :y => data[:Y],
            :a => data[:a0],
            :b => data[:b0],
            :c => data[:c0],
        )
        for k in 1:chains
    ]
    return model, data, inits
end

srand(12345)
dist_true = MixtureModel(Normal[Normal(-1.0, 1.0), Normal(0.0,1.0), Normal(1.0,1.0)], [0.2, 0.4, 0.4])
n = 2^7
Y = rand(dist_true, n)
beta = 1.0

println("========== Sample")
@show dist_true
@show n
@show mean(Y)
@show std(Y)
println("========== beta")
@show beta
println("="^20)

iterations = 4200
burnin = 200
thin = 1
chains = 8
model, data, inits = sample2model_mixnormal1(Y, beta=beta, chains=chains)
@time sim = mcmc(model, data, inits, iterations, burnin=burnin, thin=thin, chains=chains, verbose=false)
showsummary(sim, show_describe=true, show_gelman=false, plot_trace=true, plot_chain=true)

通常の正規分布モデルnormal¶

In [ ]:

@everywhere normal(mu, sigma) = Normal(mu, sigma)
@everywhere thermal_normal(mu, sigma, beta) = ThermalCU(normal(mu, sigma), beta)

normal(w::AbstractVector) = normal(w[1], w[2])
unlink_normal(w) = [w[1], log(w[2])]  # unlink_normal : R×(0,∞) -> R^2
link_normal(z)   = [z[2], exp(z[2])]  #   link_normal : R^2 → R×(0,∞)

pdf_normal(w,x)  = pdf(normal(w),x)
function rand_prior_normal(M)
    local mu = randn(M)
    local sigma = abs.(randn(M))
    return [mu';sigma']
end

function sample2model_normal(Y;
        beta = 1.0,
        dist_model = thermal_normal,
        mu0    = 0.0,
        sigma0 = 1.0,
        prior_mu = Normal(0.0, 1.0),
        prior_sigma = Truncated(Normal(1.0, 1.0), 0, Inf),
        chains = 2
    )
    local data = Dict{Symbol, Any}(
        :Y => Y,
        :n => length(Y),
        :mu0 => mu0,
        :sigma0 => sigma0,
        :beta => beta,
    )
    
    local model = Model(
        y = Stochastic(1, (mu, sigma) -> dist_model(mu, sigma, data[:beta]), false),
        mu = Stochastic(() -> prior_mu, true),
        sigma = Stochastic(() -> prior_sigma, true),
    )
    local scheme = [
        NUTS([:mu, :sigma])
    ]
    setsamplers!(model, scheme)
    
    local inits = [
        Dict{Symbol, Any}(
            :y => data[:Y],
            :mu => data[:mu0],
            :sigma => data[:sigma0],
        )
        for k in 1:chains
    ]
    return model, data, inits
end

srand(12345)
dist_true = Normal(0.5,1.0)
n = 2^7
Y = rand(dist_true, n)
beta = 1.0

println("========== Sample")
@show dist_true
@show n
@show mean(Y)
@show std(Y)
println("========== beta")
@show beta
println("="^20)

iterations = 4200
burnin = 200
thin = 1
chains = 4
model, data, inits = sample2model_normal(Y, beta=beta, chains=chains)
@time sim_normal_1 = mcmc(model, data, inits, iterations, burnin=burnin, thin=thin, chains=chains, verbose=false)
sim = sim_normal_1
showsummary(sim, show_describe=true, show_gelman=false, plot_trace=true, plot_chain=true)

In [ ]:

@everywhere normal(mu, sigma) = Normal(mu, sigma)
@everywhere thermal_normal(mu, sigma, beta) = ThermalCU(normal(mu, sigma), beta)

normal(w::AbstractVector) = normal(w[1], w[2])
unlink_normal(w) = [w[1], log(w[2])]  # unlink_normal : R×(0,∞) -> R^2
link_normal(z)   = [z[2], exp(z[2])]  #   link_normal : R^2 → R×(0,∞)

pdf_normal(w,x)  = pdf(normal(w),x)
function rand_prior_normal(M)
    local mu = randn(M)
    local sigma = abs.(randn(M))
    return [mu';sigma']
end

function sample2model_normal(Y;
        beta = 1.0,
        dist_model = thermal_normal,
        mu0    = 0.0,
        sigma0 = 1.0,
        prior_mu = Normal(0.0, 1.0),
        prior_sigma = Truncated(Normal(1.0, 1.0), 0, Inf),
        chains = 2
    )
    local data = Dict{Symbol, Any}(
        :Y => Y,
        :n => length(Y),
        :mu0 => mu0,
        :sigma0 => sigma0,
        :beta => beta,
    )
    
    local model = Model(
        y = Stochastic(1, (mu, sigma) -> dist_model(mu, sigma, data[:beta]), false),
        mu = Stochastic(() -> prior_mu, true),
        sigma = Stochastic(() -> prior_sigma, true),
    )
    local scheme = [
        NUTS([:mu, :sigma])
    ]
    setsamplers!(model, scheme)
    
    local inits = [
        Dict{Symbol, Any}(
            :y => data[:Y],
            :mu => data[:mu0],
            :sigma => data[:sigma0],
        )
        for k in 1:chains
    ]
    return model, data, inits
end

srand(12345)
dist_true = Normal(0.5,1.0)
n = 2^7
Y = rand(dist_true, n)
beta = 50.0

println("========== Sample")
@show dist_true
@show n
@show mean(Y)
@show std(Y)
println("========== beta")
@show beta
println("="^20)

iterations = 4200
burnin = 200
thin = 1
chains = 4
model, data, inits = sample2model_normal(Y, beta=beta, chains=chains)
@time sim_normal_50 = mcmc(model, data, inits, iterations, burnin=burnin, thin=thin, chains=chains, verbose=false)
sim = sim_normal_50
showsummary(sim, show_describe=true, show_gelman=false, plot_trace=true, plot_chain=true)

以下は雑多な試み¶

In [ ]:

@everywhere mutable struct ThermalMixNormal{S<:Real, T<:AbstractArray{S,1}} <:ContinuousUnivariateDistribution
    a::T
    b::T
    β::S
end

function Distributions.pdf(d::ThermalMixNormal, x::Real)
    local p = sum(d.a .* exp.(-0.5.*(x .- d.b).^2))/sqrt(2*pi)
    return p^d.β
end

@everywhere function Distributions.logpdf(d::ThermalMixNormal, x::Real)
    local l = log(sum(d.a .* exp.(-0.5.*(x .- d.b).^2))) - 0.5*log(2*pi)
    return d.β*l
end

@everywhere Distributions.maximum(d::ThermalMixNormal) = Inf
@everywhere Distributions.minimum(d::ThermalMixNormal) = -Inf

function Distributions.rand(d::ThermalMixNormal)
    local i = rand(Categorical(d.a))
    rand(Normal(d.b[i], 1.0))
end

function sample2model_ThermalMixNormal(Y;
        beta = 1.0,
        dist_model = ThermalMixNormal,
        a0 = [0.33, 0.33, 0.34],
        b0 = zeros(3),
        prior_a = Dirichlet(3, 1.0),
        prior_b = MvNormal(zeros(3), 4*ones(3)),
        chains = 2
    )
    local data = Dict{Symbol, Any}(
        :Y => Y,
        :n => length(Y),
        :a0 => a0,
        :b0 => b0,
        :beta => beta,
    )
    
    local model = Model(
        y = Stochastic(1, (a, b) -> dist_model(a, b, data[:beta]), false),
        a = Stochastic(1, () -> prior_a, true),
        b = Stochastic(1, () -> prior_b, true),
    )
    local scheme = [
        SliceSimplex(:a),
        NUTS(:b)
    ]
    setsamplers!(model, scheme)
    
    local inits = [
        Dict{Symbol, Any}(
            :y => data[:Y],
            :a => data[:a0],
            :b => data[:b0],
        )
        for k in 1:chains
    ]
    return model, data, inits
end

srand(12345)
dist_true = ThermalMixNormal([0.2, 0.3, 0.5], [-3.0, 0.0, 3.0], 1.0)
n = 2^8
Y = rand(dist_true, n)
beta = 1.0

println("========== Sample")
@show dist_true
@show n
@show mean(Y)
@show std(Y)
println("========== beta")
@show beta
println("="^20)

iterations = 1200
burnin = 200
thin = 1
chains = 8
model, data, inits = sample2model_ThermalMixNormal(Y, beta=beta, chains=chains)
@time sim = mcmc(model, data, inits, iterations, burnin=burnin, thin=thin, chains=chains, verbose=false)
showsummary(sim, show_describe=false, show_gelman=true, plot_trace=true, plot_chain=true, maxfigsize=12)

a_all = vcat((vec(sim[:,name,:].value) for name in ["a[1]", "a[2]", "a[3]"])...)
b_all = vcat((vec(sim[:,name,:].value) for name in ["b[1]", "b[2]", "b[3]"])...)
chain = hcat(a_all, b_all)
paramname = ["a", "b"]
plotchain(chain, paramname)

function makefunc_pred_TMN(sim)
    local chain = sim2chain(sim, sort(keys(sim)))
    local L = size(chain, 1)
    local a = [chain[l,1:3] for l in 1:L]
    local b = [chain[l,4:6] for l in 1:L]
    local pred_TMN(x) = sum(pdf(ThermalMixNormal(a[l], b[l], 1.0), x) for l in 1:L)/L
    return pred_TMN
end
pred_TMNs = makefunc_pred_TMN(sim)

plt.figure(figsize=(6,4.2))
x = linspace(-7,7,201)
plt.plot(x, pdf.(dist_true, x), label="true dist")
plt.plot(x, pred_TMNs.(x), label="predictive", ls="--")
plt.grid(ls=":")
plt.legend()
plt.title("Sample size n = $n")

In [ ]:

@everywhere mutable struct ThermalMixNormal{S<:Real, T<:AbstractArray{S,1}} <:ContinuousUnivariateDistribution
    a::T
    b::T
    β::S
end

function Distributions.pdf(d::ThermalMixNormal, x::Real)
    local p = sum(d.a .* exp.(-0.5.*(x .- d.b).^2))/sqrt(2*pi)
    return p^d.β
end

@everywhere function Distributions.logpdf(d::ThermalMixNormal, x::Real)
    local l = log(sum(d.a .* exp.(-0.5.*(x .- d.b).^2))) - 0.5*log(2*pi)
    return d.β*l
end

@everywhere Distributions.maximum(d::ThermalMixNormal) = Inf
@everywhere Distributions.minimum(d::ThermalMixNormal) = -Inf

function Distributions.rand(d::ThermalMixNormal)
    local i = rand(Categorical(d.a))
    rand(Normal(d.b[i], 1.0))
end

function sample2model_ThermalMixNormal(Y;
        beta = 1.0,
        dist_model = ThermalMixNormal,
        a0 = [0.33, 0.33, 0.34],
        b0 = zeros(3),
        prior_a = Dirichlet(3, 1.0),
        prior_b = MvNormal(zeros(3), 4*ones(3)),
        chains = 2
    )
    local data = Dict{Symbol, Any}(
        :Y => Y,
        :n => length(Y),
        :a0 => a0,
        :b0 => b0,
        :beta => beta,
    )
    
    local model = Model(
        y = Stochastic(1, (a, b) -> dist_model(a, b, data[:beta]), false),
        a = Stochastic(1, () -> prior_a, true),
        b = Stochastic(1, () -> prior_b, true),
    )
    local scheme = [
        SliceSimplex(:a),
        NUTS(:b)
    ]
    setsamplers!(model, scheme)
    
    local inits = [
        Dict{Symbol, Any}(
            :y => data[:Y],
            :a => data[:a0],
            :b => data[:b0],
        )
        for k in 1:chains
    ]
    return model, data, inits
end

srand(12345)
dist_true = ThermalMixNormal([0.2, 0.3, 0.5], [-3.0, 0.0, 3.0], 1.0)
n = 2^10
Y = rand(dist_true, n)
beta = 1.0

println("========== Sample")
@show dist_true
@show n
@show mean(Y)
@show std(Y)
println("========== beta")
@show beta
println("="^20)

iterations = 1200
burnin = 200
thin = 1
chains = 8
model, data, inits = sample2model_ThermalMixNormal(Y, beta=beta, chains=chains)
@time sim = mcmc(model, data, inits, iterations, burnin=burnin, thin=thin, chains=chains, verbose=false)
showsummary(sim, show_describe=false, show_gelman=true, plot_trace=true, plot_chain=true, maxfigsize=12)

a_all = vcat((vec(sim[:,name,:].value) for name in ["a[1]", "a[2]", "a[3]"])...)
b_all = vcat((vec(sim[:,name,:].value) for name in ["b[1]", "b[2]", "b[3]"])...)
chain = hcat(a_all, b_all)
paramname = ["a", "b"]
plotchain(chain, paramname)

function makefunc_pred_TMN(sim)
    local chain = sim2chain(sim, sort(keys(sim)))
    local L = size(chain, 1)
    local a = [chain[l,1:3] for l in 1:L]
    local b = [chain[l,4:6] for l in 1:L]
    local pred_TMN(x) = sum(pdf(ThermalMixNormal(a[l], b[l], 1.0), x) for l in 1:L)/L
    return pred_TMN
end
pred_TMNs = makefunc_pred_TMN(sim)

plt.figure(figsize=(6,4.2))
x = linspace(-7,7,201)
plt.plot(x, pdf.(dist_true, x), label="true dist")
plt.plot(x, pred_TMNs.(x), label="predictive", ls="--")
plt.grid(ls=":")
plt.legend()
plt.title("Sample size n = $n")

In [ ]:

@everywhere mutable struct ThermalMixNormal{S<:Real, T<:AbstractArray{S,1}} <:ContinuousUnivariateDistribution
    a::T
    b::T
    β::S
end

function Distributions.pdf(d::ThermalMixNormal, x::Real)
    local p = sum(d.a .* exp.(-0.5.*(x .- d.b).^2))/sqrt(2*pi)
    return p^d.β
end

@everywhere function Distributions.logpdf(d::ThermalMixNormal, x::Real)
    local l = log(sum(d.a .* exp.(-0.5.*(x .- d.b).^2))) - 0.5*log(2*pi)
    return d.β*l
end

@everywhere Distributions.maximum(d::ThermalMixNormal) = Inf
@everywhere Distributions.minimum(d::ThermalMixNormal) = -Inf

function Distributions.rand(d::ThermalMixNormal)
    local i = rand(Categorical(d.a))
    rand(Normal(d.b[i], 1.0))
end

function sample2model_ThermalMixNormal(Y;
        beta = 1.0,
        dist_model = ThermalMixNormal,
        a0 = [0.33, 0.33, 0.34],
        b0 = zeros(3),
        prior_a = Dirichlet(3, 1.0),
        prior_b = MvNormal(zeros(3), 4*ones(3)),
        chains = 2
    )
    local data = Dict{Symbol, Any}(
        :Y => Y,
        :n => length(Y),
        :a0 => a0,
        :b0 => b0,
        :beta => beta,
    )
    
    local model = Model(
        y = Stochastic(1, (a, b) -> dist_model(a, b, data[:beta]), false),
        a = Stochastic(1, () -> prior_a, true),
        b = Stochastic(1, () -> prior_b, true),
    )
    local scheme = [
        SliceSimplex(:a),
        NUTS(:b)
    ]
    setsamplers!(model, scheme)
    
    local inits = [
        Dict{Symbol, Any}(
            :y => data[:Y],
            :a => data[:a0],
            :b => data[:b0],
        )
        for k in 1:chains
    ]
    return model, data, inits
end

srand(12345)
dist_true = ThermalMixNormal([0.35, 0.65], [-1.5, 1.5], 1.0)
n = 2^8
Y = rand(dist_true, n)
beta = 1.0

println("========== Sample")
@show dist_true
@show n
@show mean(Y)
@show std(Y)
println("========== beta")
@show beta
println("="^20)

iterations = 1200
burnin = 200
thin = 1
chains = 8
model, data, inits = sample2model_ThermalMixNormal(Y, beta=beta, chains=chains)
@time sim = mcmc(model, data, inits, iterations, burnin=burnin, thin=thin, chains=chains, verbose=false)
showsummary(sim, show_describe=false, show_gelman=true, plot_trace=true, plot_chain=true, maxfigsize=12)

a_all = vcat((vec(sim[:,name,:].value) for name in ["a[1]", "a[2]", "a[3]"])...)
b_all = vcat((vec(sim[:,name,:].value) for name in ["b[1]", "b[2]", "b[3]"])...)
chain = hcat(a_all, b_all)
paramname = ["a", "b"]
plotchain(chain, paramname)

function makefunc_pred_TMN(sim)
    local chain = sim2chain(sim, sort(keys(sim)))
    local L = size(chain, 1)
    local a = [chain[l,1:3] for l in 1:L]
    local b = [chain[l,4:6] for l in 1:L]
    local pred_TMN(x) = sum(pdf(ThermalMixNormal(a[l], b[l], 1.0), x) for l in 1:L)/L
    return pred_TMN
end
pred_TMNs = makefunc_pred_TMN(sim)

plt.figure(figsize=(6,4.2))
x = linspace(-7,7,201)
plt.plot(x, pdf.(dist_true, x), label="true dist")
plt.plot(x, pred_TMNs.(x), label="predictive", ls="--")
plt.grid(ls=":")
plt.legend()
plt.title("Sample size n = $n")

In [ ]:

@everywhere mutable struct ThermalMixNormal{S<:Real, T<:AbstractArray{S,1}} <:ContinuousUnivariateDistribution
    a::T
    b::T
    β::S
end

function Distributions.pdf(d::ThermalMixNormal, x::Real)
    local p = sum(d.a .* exp.(-0.5.*(x .- d.b).^2))/sqrt(2*pi)
    return p^d.β
end

@everywhere function Distributions.logpdf(d::ThermalMixNormal, x::Real)
    local l = log(sum(d.a .* exp.(-0.5.*(x .- d.b).^2))) - 0.5*log(2*pi)
    return d.β*l
end

@everywhere Distributions.maximum(d::ThermalMixNormal) = Inf
@everywhere Distributions.minimum(d::ThermalMixNormal) = -Inf

function Distributions.rand(d::ThermalMixNormal)
    local i = rand(Categorical(d.a))
    rand(Normal(d.b[i], 1.0))
end

function sample2model_ThermalMixNormal(Y;
        beta = 1.0,
        dist_model = ThermalMixNormal,
        a0 = [0.33, 0.33, 0.34],
        b0 = zeros(3),
        prior_a = Dirichlet(3, 1.0),
        prior_b = MvNormal(zeros(3), 4*ones(3)),
        chains = 2
    )
    local data = Dict{Symbol, Any}(
        :Y => Y,
        :n => length(Y),
        :a0 => a0,
        :b0 => b0,
        :beta => beta,
    )
    
    local model = Model(
        y = Stochastic(1, (a, b) -> dist_model(a, b, data[:beta]), false),
        a = Stochastic(1, () -> prior_a, true),
        b = Stochastic(1, () -> prior_b, true),
    )
    local scheme = [
        SliceSimplex(:a),
        NUTS(:b)
    ]
    setsamplers!(model, scheme)
    
    local inits = [
        Dict{Symbol, Any}(
            :y => data[:Y],
            :a => data[:a0],
            :b => data[:b0],
        )
        for k in 1:chains
    ]
    return model, data, inits
end

srand(12345)
dist_true = ThermalMixNormal([0.35, 0.65], [-1.5, 1.5], 1.0)
n = 2^10
Y = rand(dist_true, n)
beta = 1.0

println("========== Sample")
@show dist_true
@show n
@show mean(Y)
@show std(Y)
println("========== beta")
@show beta
println("="^20)

iterations = 1200
burnin = 200
thin = 1
chains = 8
model, data, inits = sample2model_ThermalMixNormal(Y, beta=beta, chains=chains)
@time sim = mcmc(model, data, inits, iterations, burnin=burnin, thin=thin, chains=chains, verbose=false)
showsummary(sim, show_describe=false, show_gelman=true, plot_trace=true, plot_chain=true, maxfigsize=12)

a_all = vcat((vec(sim[:,name,:].value) for name in ["a[1]", "a[2]", "a[3]"])...)
b_all = vcat((vec(sim[:,name,:].value) for name in ["b[1]", "b[2]", "b[3]"])...)
chain = hcat(a_all, b_all)
paramname = ["a", "b"]
plotchain(chain, paramname)

function makefunc_pred_TMN(sim)
    local chain = sim2chain(sim, sort(keys(sim)))
    local L = size(chain, 1)
    local a = [chain[l,1:3] for l in 1:L]
    local b = [chain[l,4:6] for l in 1:L]
    local pred_TMN(x) = sum(pdf(ThermalMixNormal(a[l], b[l], 1.0), x) for l in 1:L)/L
    return pred_TMN
end
pred_TMNs = makefunc_pred_TMN(sim)

plt.figure(figsize=(6,4.2))
x = linspace(-7,7,201)
plt.plot(x, pdf.(dist_true, x), label="true dist")
plt.plot(x, pred_TMNs.(x), label="predictive", ls="--")
plt.grid(ls=":")
plt.legend()
plt.title("Sample size n = $n")

In [ ]:

# 分散を固定しない混合正規分布モデル

##### ThermalMixtureNormal <:ContinuousUnivariateDistribution

@everywhere begin
    mutable struct ThermalMixtureNormal{S<:Real, T<:AbstractArray{S,1}} <: ContinuousUnivariateDistribution
        μ::T
        σ::T
        p::T
        β::S
    end
end

@everywhere begin
    import Distributions: pdf, logpdf, maximum, minimum, rand, _logpdf, length, insupport, _rand!   

    function pdf(d::ThermalMixtureNormal, x::Real)
        local s = sum( d.p .* exp.(-(x .- d.μ).^2./(2*d.σ.^2)) ./ sqrt.(2*pi*d.σ.^2) )
        return s^d.β
    end

    function logpdf(d::ThermalMixtureNormal, x::Real)
        local l = log(sum( d.p .* exp.(-(x .- d.μ).^2./(2*d.σ.^2)) ./ sqrt.(2*pi*d.σ.^2) ))
        return d.β*l
    end

    minimum(d::ThermalMixtureNormal) = -Inf
    maximum(d::ThermalMixtureNormal) =  Inf

    function rand(d::ThermalMixtureNormal)
        local i = rand(Categorical(d.p))
        rand(Normal(d.μ[i], d.σ[i]))
    end
end

##### MvHalfNormal <:ContinuousMultivariateDistribution

@everywhere begin
    mutable struct ThermalMvHalfNormal{S<:Real, T<:AbstractArray{S,1}} <: ContinuousMultivariateDistribution
        σ::T
        β::S
    end
end

@everywhere begin
    import Distributions: pdf, logpdf, maximum, minimum, rand, _logpdf, length, insupport, _rand!   

    ThermalMvHalfNormal(k::T, β::S) where {T<:Integer, S<:Real} = ThermalMvHalfNormal(ones(k), β)
    ThermalMvHalfNormal(k::T, σ::S, β::S) where {T<:Integer, S<:Real} = ThermalMvHalfNormal(σ*ones(k), β)

    function pdf(d::ThermalMvHalfNormal, x::AbstractArray{S,1}) where S<:Real
        local s = prod(2 .* pdf.(Normal.(zero(S), d.σ), x))
        return s^d.β
    end

    function _logpdf(d::ThermalMvHalfNormal, x::AbstractArray{S,1}) where S<:Real
        local l = sum(log(2) .+ logpdf.(Normal.(zero(S), d.σ), x))
        return d.β*l
    end

    length(d::ThermalMvHalfNormal) = length(d.σ)

    function insupport(d::ThermalMvHalfNormal, x::AbstractArray{S,1}) where S<:Real
        all(x .> zero(S))
    end

    function _rand!(d::ThermalMvHalfNormal, x::AbstractArray{S,1}) where S<:Real
        x = abs.(rand.(Normal.(zero(S), d.σ)))
    end

    function _rand!(d::ThermalMvHalfNormal, x::AbstractArray{S,2}) where S<:Real
        for i in 1:size(x,2)
            x[:,i] = abs.(rand.(Normal.(zero(S), d.σ)))
        end
        x
    end
end

##### MCMC

function sample2model_ThermalMixtureNormal(Y;
        beta = 1.0,
        dist_model = ThermalMixtureNormal,
        μ0 = zeros(3),
        σ0 = ones(3),
        p0 = [0.33, 0.33, 0.34],
        prior_μ = MvNormal(zeros(3), 5*ones(3)),
        prior_σ = ThermalMvHalfNormal(3, 5.0, 1.0),
        prior_p = Dirichlet(3, 1.0),
        chains = 2,
    )
    local data = Dict{Symbol, Any}(
        :Y => Y,
        :n => length(Y),
        :μ0 => μ0,
        :σ0 => σ0,
        :p0 => p0,
        :beta => beta,
    )
    
    local model = Model(
        y = Stochastic(1, (μ, σ, p) -> dist_model(μ, σ, p, data[:beta]), false),
        μ = Stochastic(1, () -> prior_μ, true),
        σ = Stochastic(1, () -> prior_σ, true),
        p  = Stochastic(1, () -> prior_p,  true),
    )
    local scheme = [
        SliceSimplex(:p),
        NUTS([:μ, :σ])
    ]
    setsamplers!(model, scheme)
    
    local inits = [
        Dict{Symbol, Any}(
            :y => data[:Y],
            :μ => data[:μ0],
            :σ => data[:σ0],
            :p  => data[:p0],
        )
        for k in 1:chains
    ]
    return model, data, inits
end

srand(12345)
dist_true = ThermalMixtureNormal([-3.0, 0.0, 6.0], [1.5, 0.5, 3.0], [0.3, 0.2, 0.5], 1.0)
n = 2^8
Y = rand(dist_true, n)
beta = 1.0

println("========== Sample")
@show dist_true
@show n
@show mean(Y)
@show std(Y)
println("========== beta")
@show beta
println("="^20)

iterations = 1200
burnin = 200
thin = 1
chains = 4
model, data, inits = sample2model_ThermalMixtureNormal(Y, beta=beta, chains=chains)
@time sim = mcmc(model, data, inits, iterations, burnin=burnin, thin=thin, chains=chains, verbose=false)
showsummary(sim, show_describe=false, show_gelman=false, plot_trace=true, plot_chain=true, maxfigsize=12)

μ_all = vcat((vec(sim[:,name,:].value) for name in ["μ[1]", "μ[2]", "μ[3]"])...)
σ_all = vcat((vec(sim[:,name,:].value) for name in ["σ[1]", "σ[2]", "σ[3]"])...)
p_all = vcat((vec(sim[:,name,:].value) for name in ["p[1]", "p[2]", "p[3]"])...)
chain = hcat(p_all, μ_all, σ_all)
paramname = ["p", "μ", "σ"]
plotchain(chain, paramname)

function makefunc_pred_TMN(sim)
    local μ = sim2chain(sim, ["μ[1]", "μ[2]", "μ[3]"])
    local σ = sim2chain(sim, ["σ[1]", "σ[2]", "σ[3]"])
    local p = sim2chain(sim, ["p[1]", "p[2]", "p[3]"])
    local L = size(p,2)
    local pred_TMN(x) = sum(pdf(ThermalMixtureNormal(μ[l,:], σ[l,:], p[l,:], 1.0), x) for l in 1:L)/L
    return pred_TMN
end
pred_TMNs = makefunc_pred_TMN(sim)

plt.figure(figsize=(6,4.2))
x = linspace(-7,15,201)
plt.plot(x, pdf.(dist_true, x), label="true dist")
plt.plot(x, pred_TMNs.(x), label="predictive", ls="--")
plt.grid(ls=":")
plt.legend()
plt.title("Sample size n = $n")

In [ ]:

# 分散を固定しない混合正規分布モデル

##### ThermalMixtureNormal <:ContinuousUnivariateDistribution

@everywhere begin
    mutable struct ThermalMixtureNormal{S<:Real, T<:AbstractArray{S,1}} <: ContinuousUnivariateDistribution
        μ::T
        σ::T
        p::T
        β::S
    end
end

@everywhere begin
    import Distributions: pdf, logpdf, maximum, minimum, rand, _logpdf, length, insupport, _rand!   

    function pdf(d::ThermalMixtureNormal, x::Real)
        local s = sum( d.p .* exp.(-(x .- d.μ).^2./(2*d.σ.^2)) ./ sqrt.(2*pi*d.σ.^2) )
        return s^d.β
    end

    function logpdf(d::ThermalMixtureNormal, x::Real)
        local l = log(sum( d.p .* exp.(-(x .- d.μ).^2./(2*d.σ.^2)) ./ sqrt.(2*pi*d.σ.^2) ))
        return d.β*l
    end

    minimum(d::ThermalMixtureNormal) = -Inf
    maximum(d::ThermalMixtureNormal) =  Inf

    function rand(d::ThermalMixtureNormal)
        local i = rand(Categorical(d.p))
        rand(Normal(d.μ[i], d.σ[i]))
    end
end

##### MvHalfNormal <:ContinuousMultivariateDistribution

@everywhere begin
    mutable struct ThermalMvHalfNormal{S<:Real, T<:AbstractArray{S,1}} <: ContinuousMultivariateDistribution
        σ::T
        β::S
    end
end

@everywhere begin
    import Distributions: pdf, logpdf, maximum, minimum, rand, _logpdf, length, insupport, _rand!   

    ThermalMvHalfNormal(k::T, β::S) where {T<:Integer, S<:Real} = ThermalMvHalfNormal(ones(k), β)
    ThermalMvHalfNormal(k::T, σ::S, β::S) where {T<:Integer, S<:Real} = ThermalMvHalfNormal(σ*ones(k), β)

    function pdf(d::ThermalMvHalfNormal, x::AbstractArray{S,1}) where S<:Real
        local s = prod(2 .* pdf.(Normal.(zero(S), d.σ), x))
        return s^d.β
    end

    function _logpdf(d::ThermalMvHalfNormal, x::AbstractArray{S,1}) where S<:Real
        local l = sum(log(2) .+ logpdf.(Normal.(zero(S), d.σ), x))
        return d.β*l
    end

    length(d::ThermalMvHalfNormal) = length(d.σ)

    function insupport(d::ThermalMvHalfNormal, x::AbstractArray{S,1}) where S<:Real
        all(x .> zero(S))
    end

    function _rand!(d::ThermalMvHalfNormal, x::AbstractArray{S,1}) where S<:Real
        x = abs.(rand.(Normal.(zero(S), d.σ)))
    end

    function _rand!(d::ThermalMvHalfNormal, x::AbstractArray{S,2}) where S<:Real
        for i in 1:size(x,2)
            x[:,i] = abs.(rand.(Normal.(zero(S), d.σ)))
        end
        x
    end
end

##### MCMC

function sample2model_ThermalMixtureNormal(Y;
        beta = 1.0,
        dist_model = ThermalMixtureNormal,
        μ0 = zeros(4),
        σ0 = ones(4),
        p0 = [0.25, 0.25, 0.25, 0.25],
        prior_μ = MvNormal(zeros(4), 5*ones(4)),
        prior_σ = ThermalMvHalfNormal(4, 5.0, 1.0),
        prior_p = Dirichlet(4, 1.0),
        chains = 2,
    )
    local data = Dict{Symbol, Any}(
        :Y => Y,
        :n => length(Y),
        :μ0 => μ0,
        :σ0 => σ0,
        :p0 => p0,
        :beta => beta,
    )
    
    local model = Model(
        y = Stochastic(1, (μ, σ, p) -> dist_model(μ, σ, p, data[:beta]), false),
        μ = Stochastic(1, () -> prior_μ, true),
        σ = Stochastic(1, () -> prior_σ, true),
        p  = Stochastic(1, () -> prior_p,  true),
    )
    local scheme = [
        SliceSimplex(:p),
        NUTS([:μ, :σ])
    ]
    setsamplers!(model, scheme)
    
    local inits = [
        Dict{Symbol, Any}(
            :y => data[:Y],
            :μ => data[:μ0],
            :σ => data[:σ0],
            :p  => data[:p0],
        )
        for k in 1:chains
    ]
    return model, data, inits
end

srand(12345)
dist_true = ThermalMixtureNormal([-3.0, 0.0, 6.0], [1.5, 0.5, 3.0], [0.3, 0.2, 0.5], 1.0)
n = 2^8
Y = rand(dist_true, n)
beta = 1.0

println("========== Sample")
@show dist_true
@show n
@show mean(Y)
@show std(Y)
println("========== beta")
@show beta
println("="^20)

iterations = 1200
burnin = 200
thin = 1
chains = 4
model, data, inits = sample2model_ThermalMixtureNormal(Y, beta=beta, chains=chains)
@time sim = mcmc(model, data, inits, iterations, burnin=burnin, thin=thin, chains=chains, verbose=false)
showsummary(sim, show_describe=false, show_gelman=false, plot_trace=true, plot_chain=true, maxfigsize=16)

μ_all = vcat((vec(sim[:,name,:].value) for name in ["μ[1]", "μ[2]", "μ[3]", "μ[4]"])...)
σ_all = vcat((vec(sim[:,name,:].value) for name in ["σ[1]", "σ[2]", "σ[3]", "σ[4]"])...)
p_all = vcat((vec(sim[:,name,:].value) for name in ["p[1]", "p[2]", "p[3]", "p[4]"])...)
chain = hcat(p_all, μ_all, σ_all)
paramname = ["p", "μ", "σ"]
plotchain(chain, paramname)

function makefunc_pred_TMN(sim)
    local μ = sim2chain(sim, ["μ[1]", "μ[2]", "μ[3]", "μ[4]"])
    local σ = sim2chain(sim, ["σ[1]", "σ[2]", "σ[3]", "σ[4]"])
    local p = sim2chain(sim, ["p[1]", "p[2]", "p[3]", "p[4]"])
    local L = size(p,2)
    local pred_TMN(x) = sum(pdf(ThermalMixtureNormal(μ[l,:], σ[l,:], p[l,:], 1.0), x) for l in 1:L)/L
    return pred_TMN
end
pred_TMNs = makefunc_pred_TMN(sim)

plt.figure(figsize=(6,4.2))
x = linspace(-7,15,201)
plt.plot(x, pdf.(dist_true, x), label="true dist")
plt.plot(x, pred_TMNs.(x), label="predictive", ls="--")
plt.grid(ls=":")
plt.legend()
plt.title("Sample size n = $n")

2つの山が重なって1つの山になっているサンプル¶

In [ ]:

# 分散を固定しない混合正規分布モデル

##### ThermalMixtureNormal <:ContinuousUnivariateDistribution

@everywhere begin
    mutable struct ThermalMixtureNormal{S<:Real, T<:AbstractArray{S,1}} <: ContinuousUnivariateDistribution
        μ::T
        σ::T
        p::T
        β::S
    end
end

@everywhere begin
    import Distributions: pdf, logpdf, maximum, minimum, rand, _logpdf, length, insupport, _rand!   

    function pdf(d::ThermalMixtureNormal, x::Real)
        local s = sum( d.p .* exp.(-(x .- d.μ).^2./(2*d.σ.^2)) ./ sqrt.(2*pi*d.σ.^2) )
        return s^d.β
    end

    function logpdf(d::ThermalMixtureNormal, x::Real)
        local l = log(sum( d.p .* exp.(-(x .- d.μ).^2./(2*d.σ.^2)) ./ sqrt.(2*pi*d.σ.^2) ))
        return d.β*l
    end

    minimum(d::ThermalMixtureNormal) = -Inf
    maximum(d::ThermalMixtureNormal) =  Inf

    function rand(d::ThermalMixtureNormal)
        local i = rand(Categorical(d.p))
        rand(Normal(d.μ[i], d.σ[i]))
    end
end

##### MvHalfNormal <:ContinuousMultivariateDistribution

@everywhere begin
    mutable struct ThermalMvHalfNormal{S<:Real, T<:AbstractArray{S,1}} <: ContinuousMultivariateDistribution
        σ::T
        β::S
    end
end

@everywhere begin
    import Distributions: pdf, logpdf, maximum, minimum, rand, _logpdf, length, insupport, _rand!   

    ThermalMvHalfNormal(k::T, β::S) where {T<:Integer, S<:Real} = ThermalMvHalfNormal(ones(k), β)
    ThermalMvHalfNormal(k::T, σ::S, β::S) where {T<:Integer, S<:Real} = ThermalMvHalfNormal(σ*ones(k), β)

    function pdf(d::ThermalMvHalfNormal, x::AbstractArray{S,1}) where S<:Real
        local s = prod(2 .* pdf.(Normal.(zero(S), d.σ), x))
        return s^d.β
    end

    function _logpdf(d::ThermalMvHalfNormal, x::AbstractArray{S,1}) where S<:Real
        local l = sum(log(2) .+ logpdf.(Normal.(zero(S), d.σ), x))
        return d.β*l
    end

    length(d::ThermalMvHalfNormal) = length(d.σ)

    function insupport(d::ThermalMvHalfNormal, x::AbstractArray{S,1}) where S<:Real
        all(x .> zero(S))
    end

    function _rand!(d::ThermalMvHalfNormal, x::AbstractArray{S,1}) where S<:Real
        x = abs.(rand.(Normal.(zero(S), d.σ)))
    end

    function _rand!(d::ThermalMvHalfNormal, x::AbstractArray{S,2}) where S<:Real
        for i in 1:size(x,2)
            x[:,i] = abs.(rand.(Normal.(zero(S), d.σ)))
        end
        x
    end
end

##### MCMC

function sample2model_ThermalMixtureNormal(Y;
        beta = 1.0,
        dist_model = ThermalMixtureNormal,
        μ0 = zeros(4),
        σ0 = ones(4),
        p0 = [0.25, 0.25, 0.25, 0.25],
        prior_μ = MvNormal(zeros(4), ones(4)),
        prior_σ = ThermalMvHalfNormal(4, 1.0, 1.0),
        prior_p = Dirichlet(4, 1.0),
        chains = 2,
    )
    local data = Dict{Symbol, Any}(
        :Y => Y,
        :n => length(Y),
        :μ0 => μ0,
        :σ0 => σ0,
        :p0 => p0,
        :beta => beta,
    )
    
    local model = Model(
        y = Stochastic(1, (μ, σ, p) -> dist_model(μ, σ, p, data[:beta]), false),
        μ = Stochastic(1, () -> prior_μ, true),
        σ = Stochastic(1, () -> prior_σ, true),
        p  = Stochastic(1, () -> prior_p,  true),
    )
    local scheme = [
        SliceSimplex(:p),
        NUTS([:μ, :σ])
    ]
    setsamplers!(model, scheme)
    
    local inits = [
        Dict{Symbol, Any}(
            :y => data[:Y],
            :μ => data[:μ0],
            :σ => data[:σ0],
            :p  => data[:p0],
        )
        for k in 1:chains
    ]
    return model, data, inits
end

srand(20171125)
dist_true = ThermalMixtureNormal([-1.0, 1.0], [0.8, 1.2], [0.3, 0.7], 1.0)
n = 2^8
Y = rand(dist_true, n)
beta = 1.0

println("========== Sample")
@show dist_true
@show n
@show mean(Y)
@show std(Y)
println("========== beta")
@show beta
println("="^20)

iterations = 2200
burnin = 200
thin = 1
chains = 4
model, data, inits = sample2model_ThermalMixtureNormal(Y, beta=beta, chains=chains)
@time sim = mcmc(model, data, inits, iterations, burnin=burnin, thin=thin, chains=chains, verbose=false)
showsummary(sim, show_describe=false, show_gelman=false, plot_trace=true, plot_chain=true, maxfigsize=16)

μ_all = vcat((vec(sim[:,name,:].value) for name in ["μ[1]", "μ[2]", "μ[3]", "μ[4]"])...)
σ_all = vcat((vec(sim[:,name,:].value) for name in ["σ[1]", "σ[2]", "σ[3]", "σ[4]"])...)
p_all = vcat((vec(sim[:,name,:].value) for name in ["p[1]", "p[2]", "p[3]", "p[4]"])...)
chain = hcat(p_all, μ_all, σ_all)
paramname = ["p", "μ", "σ"]
plotchain(chain, paramname)

function makefunc_pred_TMN(sim)
    local μ = sim2chain(sim, ["μ[1]", "μ[2]", "μ[3]", "μ[4]"])
    local σ = sim2chain(sim, ["σ[1]", "σ[2]", "σ[3]", "σ[4]"])
    local p = sim2chain(sim, ["p[1]", "p[2]", "p[3]", "p[4]"])
    local L = size(p,2)
    local pred_TMN(x) = sum(pdf(ThermalMixtureNormal(μ[l,:], σ[l,:], p[l,:], 1.0), x) for l in 1:L)/L
    return pred_TMN
end
pred_TMN = makefunc_pred_TMN(sim)
kde_sample = makefunc_pdfkde(Y)

plt.figure(figsize=(6,4.2))
x = linspace(-5,5,201)
plt.plot(x, pdf.(dist_true, x), label="true dist")
plt.plot(x, pred_TMN.(x), label="predictive", ls="--")
plt.plot(x, kde_sample.(x), label="sample kde", ls="-.")
plt.grid(ls=":")
plt.legend()
plt.title("Sample size n = $n")

In [ ]:

# 分散を固定しない混合正規分布モデル

##### ThermalMixtureNormal <:ContinuousUnivariateDistribution

@everywhere begin
    mutable struct ThermalMixtureNormal{S<:Real, T<:AbstractArray{S,1}} <: ContinuousUnivariateDistribution
        μ::T
        σ::T
        p::T
        β::S
    end
end

@everywhere begin
    import Distributions: pdf, logpdf, maximum, minimum, rand, _logpdf, length, insupport, _rand!   

    function pdf(d::ThermalMixtureNormal, x::Real)
        local s = sum( d.p .* exp.(-(x .- d.μ).^2./(2*d.σ.^2)) ./ sqrt.(2*pi*d.σ.^2) )
        return s^d.β
    end

    function logpdf(d::ThermalMixtureNormal, x::Real)
        local l = log(sum( d.p .* exp.(-(x .- d.μ).^2./(2*d.σ.^2)) ./ sqrt.(2*pi*d.σ.^2) ))
        return d.β*l
    end

    minimum(d::ThermalMixtureNormal) = -Inf
    maximum(d::ThermalMixtureNormal) =  Inf

    function rand(d::ThermalMixtureNormal)
        local i = rand(Categorical(d.p))
        rand(Normal(d.μ[i], d.σ[i]))
    end
end

##### MvHalfNormal <:ContinuousMultivariateDistribution

@everywhere begin
    mutable struct ThermalMvHalfNormal{S<:Real, T<:AbstractArray{S,1}} <: ContinuousMultivariateDistribution
        σ::T
        β::S
    end
end

@everywhere begin
    import Distributions: pdf, logpdf, maximum, minimum, rand, _logpdf, length, insupport, _rand!   

    ThermalMvHalfNormal(k::T, β::S) where {T<:Integer, S<:Real} = ThermalMvHalfNormal(ones(k), β)
    ThermalMvHalfNormal(k::T, σ::S, β::S) where {T<:Integer, S<:Real} = ThermalMvHalfNormal(σ*ones(k), β)

    function pdf(d::ThermalMvHalfNormal, x::AbstractArray{S,1}) where S<:Real
        local s = prod(2 .* pdf.(Normal.(zero(S), d.σ), x))
        return s^d.β
    end

    function _logpdf(d::ThermalMvHalfNormal, x::AbstractArray{S,1}) where S<:Real
        local l = sum(log(2) .+ logpdf.(Normal.(zero(S), d.σ), x))
        return d.β*l
    end

    length(d::ThermalMvHalfNormal) = length(d.σ)

    function insupport(d::ThermalMvHalfNormal, x::AbstractArray{S,1}) where S<:Real
        all(x .> zero(S))
    end

    function _rand!(d::ThermalMvHalfNormal, x::AbstractArray{S,1}) where S<:Real
        x = abs.(rand.(Normal.(zero(S), d.σ)))
    end

    function _rand!(d::ThermalMvHalfNormal, x::AbstractArray{S,2}) where S<:Real
        for i in 1:size(x,2)
            x[:,i] = abs.(rand.(Normal.(zero(S), d.σ)))
        end
        x
    end
end

##### MCMC

function sample2model_ThermalMixtureNormal(Y;
        beta = 1.0,
        dist_model = ThermalMixtureNormal,
        μ0 = zeros(3),
        σ0 = ones(3),
        p0 = [0.33, 0.33, 0.34],
        prior_μ = MvNormal(zeros(3), ones(3)),
        prior_σ = ThermalMvHalfNormal(3, 1.0, 1.0),
        prior_p = Dirichlet(3, 1.0),
        chains = 2,
    )
    local data = Dict{Symbol, Any}(
        :Y => Y,
        :n => length(Y),
        :μ0 => μ0,
        :σ0 => σ0,
        :p0 => p0,
        :beta => beta,
    )
    
    local model = Model(
        y = Stochastic(1, (μ, σ, p) -> dist_model(μ, σ, p, data[:beta]), false),
        μ = Stochastic(1, () -> prior_μ, true),
        σ = Stochastic(1, () -> prior_σ, true),
        p  = Stochastic(1, () -> prior_p,  true),
    )
    local scheme = [
        SliceSimplex(:p),
        NUTS([:μ, :σ])
    ]
    setsamplers!(model, scheme)
    
    local inits = [
        Dict{Symbol, Any}(
            :y => data[:Y],
            :μ => data[:μ0],
            :σ => data[:σ0],
            :p  => data[:p0],
        )
        for k in 1:chains
    ]
    return model, data, inits
end

srand(20171125)
dist_true = ThermalMixtureNormal([-1.0, 1.0], [0.8, 1.2], [0.3, 0.7], 1.0)
n = 2^8
Y = rand(dist_true, n)
beta = 1.0

println("========== Sample")
@show dist_true
@show n
@show mean(Y)
@show std(Y)
println("========== beta")
@show beta
println("="^20)

iterations = 2200
burnin = 200
thin = 1
chains = 4
model, data, inits = sample2model_ThermalMixtureNormal(Y, beta=beta, chains=chains)
@time sim = mcmc(model, data, inits, iterations, burnin=burnin, thin=thin, chains=chains, verbose=false)
showsummary(sim, show_describe=false, show_gelman=false, plot_trace=true, plot_chain=true, maxfigsize=16)

μ_all = vcat((vec(sim[:,name,:].value) for name in ["μ[1]", "μ[2]", "μ[3]"])...)
σ_all = vcat((vec(sim[:,name,:].value) for name in ["σ[1]", "σ[2]", "σ[3]"])...)
p_all = vcat((vec(sim[:,name,:].value) for name in ["p[1]", "p[2]", "p[3]"])...)
chain = hcat(p_all, μ_all, σ_all)
paramname = ["p", "μ", "σ"]
plotchain(chain, paramname)

function makefunc_pred_TMN(sim)
    local μ = sim2chain(sim, ["μ[1]", "μ[2]", "μ[3]"])
    local σ = sim2chain(sim, ["σ[1]", "σ[2]", "σ[3]"])
    local p = sim2chain(sim, ["p[1]", "p[2]", "p[3]"])
    local L = size(p,2)
    local pred_TMN(x) = sum(pdf(ThermalMixtureNormal(μ[l,:], σ[l,:], p[l,:], 1.0), x) for l in 1:L)/L
    return pred_TMN
end
pred_TMN = makefunc_pred_TMN(sim)
kde_sample = makefunc_pdfkde(Y)

plt.figure(figsize=(6,4.2))
x = linspace(-5,5,201)
plt.plot(x, pdf.(dist_true, x), label="true dist")
plt.plot(x, pred_TMN.(x), label="predictive", ls="--")
plt.plot(x, kde_sample.(x), label="sample kde", ls="-.")
plt.grid(ls=":")
plt.legend()
plt.title("Sample size n = $n")

In [ ]:

# 分散を固定しない混合正規分布モデル

##### ThermalMixtureNormal <:ContinuousUnivariateDistribution

@everywhere begin
    mutable struct ThermalMixtureNormal{S<:Real, T<:AbstractArray{S,1}} <: ContinuousUnivariateDistribution
        μ::T
        σ::T
        p::T
        β::S
    end
end

@everywhere begin
    import Distributions: pdf, logpdf, maximum, minimum, rand, _logpdf, length, insupport, _rand!   

    function pdf(d::ThermalMixtureNormal, x::Real)
        local s = sum( d.p .* exp.(-(x .- d.μ).^2./(2*d.σ.^2)) ./ sqrt.(2*pi*d.σ.^2) )
        return s^d.β
    end

    function logpdf(d::ThermalMixtureNormal, x::Real)
        local l = log(sum( d.p .* exp.(-(x .- d.μ).^2./(2*d.σ.^2)) ./ sqrt.(2*pi*d.σ.^2) ))
        return d.β*l
    end

    minimum(d::ThermalMixtureNormal) = -Inf
    maximum(d::ThermalMixtureNormal) =  Inf

    function rand(d::ThermalMixtureNormal)
        local i = rand(Categorical(d.p))
        rand(Normal(d.μ[i], d.σ[i]))
    end
end

##### MvHalfNormal <:ContinuousMultivariateDistribution

@everywhere begin
    mutable struct ThermalMvHalfNormal{S<:Real, T<:AbstractArray{S,1}} <: ContinuousMultivariateDistribution
        σ::T
        β::S
    end
end

@everywhere begin
    import Distributions: pdf, logpdf, maximum, minimum, rand, _logpdf, length, insupport, _rand!   

    ThermalMvHalfNormal(k::T, β::S) where {T<:Integer, S<:Real} = ThermalMvHalfNormal(ones(k), β)
    ThermalMvHalfNormal(k::T, σ::S, β::S) where {T<:Integer, S<:Real} = ThermalMvHalfNormal(σ*ones(k), β)

    function pdf(d::ThermalMvHalfNormal, x::AbstractArray{S,1}) where S<:Real
        local s = prod(2 .* pdf.(Normal.(zero(S), d.σ), x))
        return s^d.β
    end

    function _logpdf(d::ThermalMvHalfNormal, x::AbstractArray{S,1}) where S<:Real
        local l = sum(log(2) .+ logpdf.(Normal.(zero(S), d.σ), x))
        return d.β*l
    end

    length(d::ThermalMvHalfNormal) = length(d.σ)

    function insupport(d::ThermalMvHalfNormal, x::AbstractArray{S,1}) where S<:Real
        all(x .> zero(S))
    end

    function _rand!(d::ThermalMvHalfNormal, x::AbstractArray{S,1}) where S<:Real
        x = abs.(rand.(Normal.(zero(S), d.σ)))
    end

    function _rand!(d::ThermalMvHalfNormal, x::AbstractArray{S,2}) where S<:Real
        for i in 1:size(x,2)
            x[:,i] = abs.(rand.(Normal.(zero(S), d.σ)))
        end
        x
    end
end

##### MCMC

function sample2model_ThermalMixtureNormal(Y;
        beta = 1.0,
        dist_model = ThermalMixtureNormal,
        μ0 = zeros(2),
        σ0 = ones(2),
        p0 = [0.5, 0.5],
        prior_μ = MvNormal(zeros(2), ones(2)),
        prior_σ = ThermalMvHalfNormal(2, 1.0, 1.0),
        prior_p = Dirichlet(2, 1.0),
        chains = 2,
    )
    local data = Dict{Symbol, Any}(
        :Y => Y,
        :n => length(Y),
        :μ0 => μ0,
        :σ0 => σ0,
        :p0 => p0,
        :beta => beta,
    )
    
    local model = Model(
        y = Stochastic(1, (μ, σ, p) -> dist_model(μ, σ, p, data[:beta]), false),
        μ = Stochastic(1, () -> prior_μ, true),
        σ = Stochastic(1, () -> prior_σ, true),
        p  = Stochastic(1, () -> prior_p,  true),
    )
    local scheme = [
        SliceSimplex(:p),
        NUTS([:μ, :σ])
    ]
    setsamplers!(model, scheme)
    
    local inits = [
        Dict{Symbol, Any}(
            :y => data[:Y],
            :μ => data[:μ0],
            :σ => data[:σ0],
            :p  => data[:p0],
        )
        for k in 1:chains
    ]
    return model, data, inits
end

srand(20171125)
dist_true = ThermalMixtureNormal([-1.0, 1.0], [0.8, 1.2], [0.3, 0.7], 1.0)
n = 2^8
Y = rand(dist_true, n)
beta = 1.0

println("========== Sample")
@show dist_true
@show n
@show mean(Y)
@show std(Y)
println("========== beta")
@show beta
println("="^20)

iterations = 2200
burnin = 200
thin = 1
chains = 4
model, data, inits = sample2model_ThermalMixtureNormal(Y, beta=beta, chains=chains)
@time sim = mcmc(model, data, inits, iterations, burnin=burnin, thin=thin, chains=chains, verbose=false)
showsummary(sim, show_describe=false, show_gelman=false, plot_trace=true, plot_chain=true, maxfigsize=16)

μ_all = vcat((vec(sim[:,name,:].value) for name in ["μ[1]", "μ[2]"])...)
σ_all = vcat((vec(sim[:,name,:].value) for name in ["σ[1]", "σ[2]"])...)
p_all = vcat((vec(sim[:,name,:].value) for name in ["p[1]", "p[2]"])...)
chain = hcat(p_all, μ_all, σ_all)
paramname = ["p", "μ", "σ"]
plotchain(chain, paramname)

function makefunc_pred_TMN(sim)
    local μ = sim2chain(sim, ["μ[1]", "μ[2]"])
    local σ = sim2chain(sim, ["σ[1]", "σ[2]"])
    local p = sim2chain(sim, ["p[1]", "p[2]"])
    local L = size(p,2)
    local pred_TMN(x) = sum(pdf(ThermalMixtureNormal(μ[l,:], σ[l,:], p[l,:], 1.0), x) for l in 1:L)/L
    return pred_TMN
end
pred_TMN = makefunc_pred_TMN(sim)
kde_sample = makefunc_pdfkde(Y)

plt.figure(figsize=(6,4.2))
x = linspace(-5,5,201)
plt.plot(x, pdf.(dist_true, x), label="true dist")
plt.plot(x, pred_TMN.(x), label="predictive", ls="--")
plt.plot(x, kde_sample.(x), label="sample kde", ls="-.")
plt.grid(ls=":")
plt.legend()
plt.title("Sample size n = $n")

In [ ]: