Engee
Личный кабинет

Расширенное использование

Определение настраиваемого распределения

Turing.jl поддерживает использование распределений из пакета Distributions.jl. Путем расширения также поддерживается использование настраиваемых распределений, которые определяются как подтипы типа Distribution из пакета Distributions.jl, а также соответствующих функций.

Ниже показан рабочий процесс определения настраиваемого распределения на примере собственной реализации простого распределения Uniform.

1. Определение типа распределения

Сначала определяется тип распределения как подтип соответствующего типа из пакета Distributions.jl.

struct CustomUniform <: ContinuousUnivariateDistribution end

2. Реализация выборки и вычисления логарифма дифференциальной функции распределения (PDF)

Затем определяются функции rand и logpdf, которые будут использоваться для выполнения модели.

# выборка из диапазона [0, 1]
Distributions.rand(rng::AbstractRNG, d::CustomUniform) = rand(rng)

# p(x) = 1 → logp(x) = 0
Distributions.logpdf(d::CustomUniform, x::Real) = zero(x)

3. Определение вспомогательных функций

Иногда может потребоваться определить ряд вспомогательных функций.

3.1. Преобразование области значений

Некоторые сэмплеры, такие как HMC, требуют неограниченной области значений априорных распределений. Поэтому для использования CustomUniform в качестве априорного распределения в модели необходимо определить также способ преобразования выборок из диапазона [0, 1] в ℝ. Для этого достаточно определить соответствующий биектор Bijector из Bijectors.jl, который Turing.jl применяет на внутреннем уровне для работы с ограниченными распределениями.

Для преобразования из диапазона [0, 1] в ℝ можно использовать биектор Logit:

Bijectors.bijector(d::CustomUniform) = Logit(0., 1.)

Для ContinuousMultivariateDistribution и ContinuousMatrixDistribution необходимо было бы сделать точно то же самое. Например, Wishart определяет распределение по определенно-положительным матрицам, поэтому bijector возвращает PDBijector при вызове с распределением Wishart в качестве аргумента. Для дискретных распределений определять биектор не нужно; по умолчанию используется биектор Identity.

Кроме того, для UnivariateDistribution можно определить minimum и maximum распределения,

Distributions.minimum(d::CustomUniform) = 0.
Distributions.maximum(d::CustomUniform) = 1.

и Bijectors.jl вернет объект Bijector по умолчанию, который называется TruncatedBijector. Он использует значения minimum и maximum для вывода правильного преобразования.

Когда необходимо преобразовать ограниченное распределение в неограниченное, например при выборке с помощью HMC, в Turing происходит следующее:

b = bijector(dist)
transformed_dist = transformed(dist, b) # результатом является распределение с преобразованной поддержкой и коррекцией для logpdf

После этого можно вызывать rand и logpdf обычным образом, где

  • rand(transformed_dist) возвращает выборку из неограниченного пространства, а

  • logpdf(transformed_dist, y) возвращает логарифмическую плотность исходного распределения, но с y в неограниченном пространстве.

Дополнительные сведения о Bijectors.jl см. в файле сведений проекта.

Обновление накопленной логарифмической вероятности в определении модели

Turing накапливает логарифмические вероятности во внутренней структуре данных, которая доступна посредством внутренней переменной __varinfo__ в определении модели (дополнительные сведения о внутренних объектах модели см. ниже). Однако поскольку пользователи не должны напрямую работать с внутренними структурами данных, им предоставляется макрос Turing.@addlogprob!, увеличивающий накопленную логарифмическую вероятность. Например, он позволяет включать в вероятность произвольные члены

using Turing

myloglikelihood(x, μ) = loglikelihood(Normal(μ, 1), x)

@model function demo(x)
    μ ~ Normal()
    Turing.@addlogprob! myloglikelihood(x, μ)
end

и отклонять выборки:

using Turing
using LinearAlgebra

@model function demo(x)
    m ~ MvNormal(zero(x), I)
    if dot(m, x) < 0
        Turing.@addlogprob! -Inf
        # Досрочное завершение вычисления модели
        return
    end

    x ~ MvNormal(m, I)
    return
end

Обратите внимание, что @addlogprob! всегда увеличивает накопленную логарифмическую вероятность независимо от предоставленного контекста выборки. Например, если вы хотите применять Turing.@addlogprob! только при вычислении логарифмической вероятности и совместной логарифмической вероятности, но не при вычислении априорного распределения модели, следует проверить тип внутренней переменной _context, например:

if DynamicPPL.leafcontext(__context__) !== Turing.PriorContext()
    Turing.@addlogprob! myloglikelihood(x, μ)
end

Внутренние объекты модели

Макрос @model принимает определение модели и преобразует его так, что вызов функции создает структуру Model для использования сэмплером. Модели можно создавать вручную без использования макроса. Если взять для примера модель gdemo, определение на основе макроса

using Turing

@model function gdemo(x)
  # Задаем априорные распределения.
  s² ~ InverseGamma(2, 3)
  m ~ Normal(0, sqrt(s²))

  # Наблюдаем каждое значение x.
  @. x ~ Normal(m, sqrt(s²))
end

model = gdemo([1.5, 2.0])

можно также реализовать (в несколько менее универсальной форме) в версии без макроса:

using Turing

# Создаем функцию модели.
function gdemo(model, varinfo, context, x)
    # Предполагаем, что s² имеет распределение InverseGamma.
    s², varinfo = DynamicPPL.tilde_assume!!(
        context,
        InverseGamma(2, 3),
        Turing.@varname(s²),
        varinfo,
    )

    # Предполагаем, что m имеет нормальное распределение.
    m, varinfo = DynamicPPL.tilde_assume!!(
        context,
        Normal(0, sqrt(s²)),
        Turing.@varname(m),
        varinfo,
    )

    # Наблюдаем каждое значение x[i] в соответствии с нормальным распределением.
    DynamicPPL.dot_tilde_observe!!(context, Normal(m, sqrt(s²)), x, Turing.@varname(x), varinfo)
end
gdemo(x) = Turing.Model(gdemo, (; x))

# Создаем экземпляр объекта Model с переменными данных.
model = gdemo([1.5, 2.0])

Копирование задач

Turing копирует задачи Julia для обеспечения эффективных алгоритмов вывода, но также предоставляет альтернативную более медленную реализацию в качестве резервного варианта. Копирование задач включено по умолчанию. Для копирования задач необходимо использовать средство TapedTask из пакета Libtask.