Интерфейс выборки

# Импортируем соответствующие библиотеки.
import AbstractMCMC
using Distributions
using Random

# Определим тип сэмплера.
struct MetropolisHastings{T, D} <: AbstractMCMC.AbstractSampler
    init_θ::T
    proposal::D
end

# Конструкторы по умолчанию.
MetropolisHastings(init_θ::Real) = MetropolisHastings(init_θ, Normal(0,1))
MetropolisHastings(init_θ::Vector{<:Real}) = MetropolisHastings(init_θ, MvNormal(zero(init_θ), I))

# Определим тип модели. Хранит функцию плотности логарифма.
struct DensityModel{F<:Function} <: AbstractMCMC.AbstractModel
    ℓπ::F
end

# Создадим очень простой тип перехода, который хранит только
# выборы параметров и логарифмическую вероятность выбора.
struct Transition{T, L}
    θ::T
    lp::L
end

# Сохраним новый выбор и его логарифмическую точность.
Transition(model::DensityModel, θ) = Transition(θ, ℓπ(model, θ))

# Определим функцию first step!, которая вызывается
# в начале выборки. Возвращаем начальный параметр, используемый
# для определения сэмплера.
function AbstractMCMC.step!(
    rng::AbstractRNG,
    model::DensityModel,
    spl::MetropolisHastings,
    N::Integer,
    ::Nothing;
    kwargs...
)
    return Transition(model, spl.init_θ)
end

# Определим функцию, которая делает базовое предположение в зависимости от одномерной
# параметризации или многомерной параметризации.
propose(spl::MetropolisHastings, model::DensityModel, θ::Real) =
    Transition(model, θ + rand(spl.proposal))
propose(spl::MetropolisHastings, model::DensityModel, θ::Vector{<:Real}) =
    Transition(model, θ + rand(spl.proposal))
propose(spl::MetropolisHastings, model::DensityModel, t::Transition) =
    propose(spl, model, t.θ)

# Вычисляет вероятность `q(θ|θcond)`, используя распределение предположений `spl.proposal`.
q(spl::MetropolisHastings, θ::Real, θcond::Real) = logpdf(spl.proposal, θ - θcond)
q(spl::MetropolisHastings, θ::Vector{<:Real}, θcond::Vector{<:Real}) =
    logpdf(spl.proposal, θ - θcond)
q(spl::MetropolisHastings, t1::Transition, t2::Transition) = q(spl, t1.θ, t2.θ)

# Вычислим плотность модели при определенной параметризации.
ℓπ(model::DensityModel, θ) = model.ℓπ(θ)
ℓπ(model::DensityModel, t::Transition) = t.lp

# Определим другую функцию шага. Возвращает переход, содержащий
# либо новое предположение (если принято), либо предыдущее предположение
# (если не принято).
function AbstractMCMC.step!(
    rng::AbstractRNG,
    model::DensityModel,
    spl::MetropolisHastings,
    ::Integer,
    θ_prev::Transition;
    kwargs...
)
    # Создадим новое предположение.
    θ = propose(spl, model, θ_prev)

    # Вычислим логарифмическую вероятность принятия.
    α = ℓπ(model, θ) - ℓπ(model, θ_prev) + q(spl, θ_prev, θ) - q(spl, θ, θ_prev)

    # Решим, вернуть ли прежнее значение θ или новое.
    if log(rand(rng)) < min(α, 0.0)
        return θ
    else
        return θ_prev
    end
end

# Базовый конструктор цепочек, работающий со структурой переходов, которую мы определили.
function AbstractMCMC.bundle_samples(
    rng::AbstractRNG,
    ℓ::DensityModel,
    s::MetropolisHastings,
    N::Integer,
    ts::Vector{<:Transition},
    chain_type::Type{Any};
    param_names=missing,
    kwargs...
)
    # Преобразуем все переходы в вектор векторов.
    vals = copy(reduce(hcat,[vcat(t.θ, t.lp) for t in ts])')

    # Проверим, получены ли имена параметров.
    if ismissing(param_names)
        param_names = ["Parameter $i" for i in 1:(length(first(vals))-1)]
    end

    # Добавим поле логарифмической плотности в имена параметров.
    push!(param_names, "lp")

    # Соберем все вместе и вернем структуру цепочек.
    return Chains(vals, param_names, (internals=["lp"],))
end

# Сгенерируем набор данных из последующих данных, которые нужно оценить.
data = rand(Normal(5, 3), 30)

# Определим компоненты базовой модели.
insupport(θ) = θ[2] >= 0
dist(θ) = Normal(θ[1], θ[2])
density(θ) = insupport(θ) ? sum(logpdf.(dist(θ), data)) : -Inf

# Создадим DensityModel.
model = DensityModel(density)

# Настроим сэмплер с начальными параметрами.
spl = MetropolisHastings([0.0, 0.0])

# Выполним выборку из последующих данных.
chain = sample(model, spl, 100000; param_names=["μ", "σ"])

Object of type Chains, with data of type 100000×3×1 Array{Float64,3}

Iterations        = 1:100000
Thinning interval = 1
Chains            = 1
Samples per chain = 100000
internals         = lp
parameters        = μ, σ

2-element Array{ChainDataFrame,1}

Summary Statistics

│ Row │ parameters │ mean    │ std      │ naive_se   │ mcse       │ ess     │ r_hat   │
│     │ Symbol     │ Float64 │ Float64  │ Float64    │ Float64    │ Any     │ Any     │
├─────┼────────────┼─────────┼──────────┼────────────┼────────────┼─────────┼─────────┤
│ 1   │ μ          │ 5.33157 │ 0.854193 │ 0.0027012  │ 0.00893069 │ 8344.75 │ 1.00009 │
│ 2   │ σ          │ 4.54992 │ 0.632916 │ 0.00200146 │ 0.00534942 │ 14260.8 │ 1.00005 │

Quantiles

│ Row │ parameters │ 2.5%    │ 25.0%   │ 50.0%   │ 75.0%   │ 97.5%   │
│     │ Symbol     │ Float64 │ Float64 │ Float64 │ Float64 │ Float64 │
├─────┼────────────┼─────────┼─────────┼─────────┼─────────┼─────────┤
│ 1   │ μ          │ 3.6595  │ 4.77754 │ 5.33182 │ 5.89509 │ 6.99651 │
│ 2   │ σ          │ 3.5097  │ 4.09732 │ 4.47805 │ 4.93094 │ 5.96821 │