Документация Engee

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

Turing реализует интерфейс выборки (размещенный в AbstractMCMC), который призван обеспечить общую основу для сэмплеров Монте-Карло цепи Маркова. Интерфейс представляет несколько структур и функций, которые необходимо перегрузить, чтобы реализовать совместимый с интерфейсом сэмплер.

В этом руководстве будет показано, как реализовать интерфейс без Turing.

Обзор интерфейса

Любая реализация метода вывода, использующая интерфейс AbstractMCMC, должна реализовать подмножество следующих типов и функций.

  1. Подтип AbstractSampler, определяемый как изменяемая структура, содержащая информацию о состоянии или параметры сэмплера.

  2. Функция sample_init!, которая выполняет любую необходимую настройку (по умолчанию настройка не выполняется).

  3. Функция step!, которая возвращает переход, представляющий собой одну выборку из сэмплера.

  4. Функция transitions_init, которая возвращает контейнер для переходов, полученных от сэмплера (по умолчанию возвращается Vector{T} длиной N, где T — этот тип перехода, полученного на первом шаге, а N — это количество запрошенных образцов).

  5. Функция transitions_save!, сохраняющая переходы в контейнер (по умолчанию переход итерации i сохраняется в позиции i в векторе переходов).

  6. Функция sample_end!, которая обрабатывает любое завершение сэмплера (по умолчанию завершение не выполняется).

  7. Функция bundle_samples, которая принимает контейнер переходов и возвращает коллекцию образцов (по умолчанию возвращается вектор переходов).

Методы интерфейса с восклицательными знаками предназначены для изменения состояния. Любая изменяемая функция призвана обеспечить изменение там, где это необходимо. Вы можете использовать следующие функции:

  • sample_init! для выполнения некоторой подготовки сэмплера перед началом выборки. Это может привести к изменению состояния сэмплера.

  • step! может изменять флаг сэмплера после каждой выборки.

  • sample_end! содержит любое завершение, которое вам может понадобиться. Если вы выполняли выборку в преобразованном пространстве, то здесь вы можете преобразовать все обратно в ограниченное пространство.

Зачем нужен интерфейс?

Мотивация использования интерфейса заключается в том, чтобы фантастическое сообщество специалистов по языку вероятностного программирования Julia могло получить набор стандартов и общих реализаций, способствующих всеобщему развитию. Методы Монте-Карло цепи Маркова, как правило, имеют очень похожую структуру, поэтому общий интерфейс должен помочь реализовать более эффективные методы вывода, встроенные в одноцелевые пакеты, в широком масштабе в рамках сообщества.

Реализация метода Метрополиса-Гастингса без Turing

Метод Метрополиса-Гастингса часто является первым методом выборки, с которым знакомятся люди. Это очень простой алгоритм и, соответственно, самый простой в реализации, поэтому он является подходящим примером. В этом разделе вы узнаете, как использовать перечисленные выше типы и функции для реализации сэмплера Метрополиса-Гастингса с помощью интерфейса MCMC.

Полный код этой реализации находится в AdvancedMH.jl.

Импорты

Начнем с импорта соответствующих библиотек. Импортируем AbstractMCMC, который содержит структуру интерфейса, которую мы будем заполнять. Нам также необходимы Distributions и Random.

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

Для расширения интерфейса (например, того, которое мы сейчас пишем) обычно требуются перегрузка или реализация нескольких функций. В частности, следует импортировать (import) функции, которые вы собираетесь перегрузить. Это выполняется в следующем блоке кода.

Из Distributions нам нужны Sampleable, VariateForm и ValueSupport — три абстрактных типа, которые определяют распределение. Предполагается, что модели в интерфейсе являются подтипами Sampleable{VariateForm, ValueSupport}. В этом разделе модель будет очень простой, поэтому мы не будем использовать их, кроме как для того, чтобы убедиться в правильности диспетчеризации функций вывода.

Сэмплер

Начнем определение сэмплера с определения сэмплера MetropolisHastings, который является подтипом AbstractSampler. Правильная типизация очень важна для надлежащей реализации интерфейса: если у вас отсутствует подтип, метод может не быть диспетчеризирован при вызове sample.

# Определим тип сэмплера.
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))

Выше мы определили сэмплер, хранящий начальную параметризацию прежнего, и объект распределения, из которого берутся предложения. У вас может быть структура, не имеющая полей. И ее можно просто использовать для диспетчеризации соответствующих функций. Или же вы можете хранить в сэмплере большое количество информации о состоянии.

Общее предположение относительно того, что хранить в структуре сэмплера, заключается в том, что в структуру сэмплера следует помещать все, что вам может понадобиться для выполнения выводов между выборками, но что вы не хотите хранить в переходе. Это единственный способ переноса информации о состоянии, не связанной с выборкой, между вызовами step!.

Модель

Далее нам понадобится некая модель. Модель — это структура, являющаяся подтипом AbstractModel, которая содержит любую информацию, необходимую для выполнения выводов для имеющейся задачи. В нашем случае требуется знать параметры среднего и дисперсии для стандартного нормального распределения, поэтому можно свести модель к логарифмической плотности нормального распределения.

Обратите внимание, что это приходится делать только потому, что мы еще не интегрировали сэмплер с Turing. Turing имеет очень сложный механизм моделирования, который избавляет от необходимости определять пользовательские структуры модели.

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

Переход

Следующим шагом является определение перехода, который мы будем возвращать из каждого вызова step!. Мы упростим задачу, просто определив структуру-оболочку, которая будет содержать выборы параметров и логарифмическую плотность этого выбора:

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

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

Transition теперь может хранить любой тип параметра, будь то вектор выборов из нескольких параметров или один одномерный выбор.

Схема Метрополиса-Гастингса

Теперь пришло время заняться выводом. Мы определили все необходимые основные части, но нам нужно реализовать функцию step!, которая на самом деле выполняет вывод.

Вспомним, что схема Метрополиса-Гастингса реализует очень простой алгоритм:

  1. Выбор начального состояния .

  2. Для в ,]

    • Генерирование параметризации предложения

    • Вычисление вероятности принятия,

    • Если , где ,], то Иначе,

Конечно, гораздо проще сделать это в логарифмическом пространстве, поэтому вероятность принятия чаще всего записывается следующим образом.

С точки зрения интерфейса нужно сделать следующее.

  1. Создать новый переход, содержащий предложенный образец.

  2. Вычислить вероятность принятия решения.

  3. Если мы принимаем решение, возвратить новый переход, в противном случае возвратить старый.

Шаги

Функция step! выполняет основную часть вывода. В нашем случае мы реализуем две функции step! — одну для самой первой итерации и одну для каждой последующей.

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

Первая функция step! просто упаковывает начальную параметризацию внутри сэмплера и возвращает ее. Мы неявным образом принимаем самую первую параметризацию.

Другая функция step! выполняет обычные шаги алгоритма Метрополиса-Гастингса. В комплект входят несколько вспомогательных функций proposal и q, которые должны реплицировать функции в приведенном выше псевдокоде.

  • proposal генерирует новое предположение в виде Transition, которое может быть одномерным, если передаваемое значение является одномерным, или многомерным, если заданное значение Transition является многомерным. Предположения используют базовое распределение Normal или MvNormal.

  • q возвращает логарифмическую плотность одной параметризации, обусловленную другой, в соответствии с распределением предположения.

  • step! генерирует новое предположение, проверяет вероятность принятия, а затем возвращает либо предыдущий переход, либо предположительный переход.

# Определим функцию, которая делает базовое предположение в зависимости от одномерной
# параметризации или многомерной параметризации.
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

Цепочки

В реализации по умолчанию sample просто возвращает вектор всех переходов. Если нужно получить объект Chains (например, для упрощения последующего анализа), также необходимо реализовать функцию bundle_samples. Она принимает вектор переходов и возвращает коллекцию образцов. К счастью, переход (Transition) невероятно прост, и нужно создать лишь небольшую функциональность для принятия пользовательских имен параметров, передаваемых пользователем.

# Базовый конструктор цепочек, работающий со структурой переходов, которую мы определили.
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

Все готово.

Вы даже можете реализовать различные форматы вывода, реализовав bundle_samples для различных chain_type, которые могут быть указаны в качестве именованного аргумента для sample. По умолчанию sample использует chain_type = Any.

Тестирование реализации

Теперь, когда у нас есть все части, необходимо протестировать реализацию, определив модель для вычисления параметров среднего и дисперсии нормального распределения. Это можно сделать, построив целевую функцию плотности, предоставив выборку данных, а затем запустив сэмплер с помощью sample.

# Сгенерируем набор данных из последующих данных, которые нужно оценить.
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=["μ", "σ"])

Если все функции интерфейса были расширены должным образом, вы должны получить вывод из display(chain), который имеет следующий вид:

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 │

Похоже, мы очень близки к истинным параметрам Normal(5,3), хотя и с довольно большой дисперсией из-за малого размера выборки.

Заключение

Мы рассмотрели, как реализовать интерфейс выборки для общих проектов. Методы интерфейса Turing постоянно развиваются, поэтому открывайте проблемы по адресу AbstractMCMC с запросами на получение функций или для решения возникающих вопросов.