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

Визуализация сэмплера

Код

Для каждого сэмплера мы будем использовать один и тот же код для построения графиков его путей. Приведенный ниже блок загружает соответствующие библиотеки и определяет функцию для построения графика траектории сэмплера по апостериорному распределению.

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

ENV["GKS_ENCODING"] = "utf-8" # Позволяет использовать символы Unicode в Plots.j.
using Plots
using StatsPlots
using Turing
using Bijectors
using Random
using DynamicPPL: getlogp, settrans!, getval, reconstruct, vectorize, setval!

# Зададим начальное значение.
Random.seed!(0)

# Определим необычную модель.
@model function gdemo(x)
    s² ~ InverseGamma(2, 3)
    m ~ Normal(0, sqrt(s²))
    bumps = sin(m) + cos(m)
    m = m + 5*bumps
    for i in eachindex(x)
      x[i] ~ Normal(m, sqrt(s²))
    end
    return s², m
end

# Определим точки данных.
x = [1.5, 2.0, 13.0, 2.1, 0.0]

# Настроим вызов модели, сделаем выборку из априорного распределения.
model = gdemo(x)
vi = Turing.VarInfo(model)

# Перед выборкой преобразуем параметр дисперсии в вещественную строку.
# Примечание. Нам приходится делать это только здесь, потому что мы выполняем практическое задание.
# Во время обычной выборки Turing сделает все это за вас.
dist = InverseGamma(2,3)
svn = vi.metadata.s².vns[1]
mvn = vi.metadata.m.vns[1]
setval!(vi, vectorize(dist, Bijectors.link(dist, reconstruct(dist, getval(vi, svn)))), svn)
settrans!(vi, true, svn)

# Оценим поверхность по координатам.
function evaluate(m1, m2)
    spl = Turing.SampleFromPrior()
    vi[svn] = [m1]
    vi[mvn] = [m2]
    model(vi, spl)
    getlogp(vi)
end

function plot_sampler(chain; label="")
    # Извлечем значения из цепочки.
    val = get(chain, [:s², :m, :lp])
    ss = link.(Ref(InverseGamma(2, 3)), val.s²)
    ms = val.m
    lps = val.lp

    # Количество точек поверхности для выборки.
    granularity = 100

    # Точки начала/остановки диапазона.
    spread = 0.5
    σ_start = minimum(ss) - spread * std(ss); σ_stop = maximum(ss) + spread * std(ss);
    μ_start = minimum(ms) - spread * std(ms); μ_stop = maximum(ms) + spread * std(ms);
    σ_rng = collect(range(σ_start, stop=σ_stop, length=granularity))
    μ_rng = collect(range(μ_start, stop=μ_stop, length=granularity))

    # Построим график поверхности.
    p = surface(σ_rng, μ_rng, evaluate,
          camera=(30, 65),
        #   ticks=nothing,
          colorbar=false,
          color=:inferno,
          title=label)

    line_range = 1:length(ms)

    scatter3d!(ss[line_range], ms[line_range], lps[line_range],
        mc =:viridis, marker_z=collect(line_range), msw=0,
        legend=false, colorbar=false, alpha=0.5,
        xlabel="σ", ylabel="μ", zlabel="Log probability",
        title=label)

    return p
end;

Сэмплеры

Гиббс

Выборка Гиббса обычно демонстрирует траекторию с отклонениями. В приведенном ниже примере показано сочетание выборки HMC и PG для прохода по апостериорному распределению.

c = sample(model, Gibbs(HMC(0.01, 5, :s), PG(20, :m)), 1000)
plot_sampler(c)
samplers 1

HMC

Выборка гамильтониана Монте-Карло (HMC) — это типичный используемый сэмплер, поскольку он достаточно хорошо сходится. Однако часто бывает сложно задать правильные параметры для этого сэмплера, и сэмплер NUTS часто проще запустить, если вы не хотите тратить много времени на определение размера шага и количества шагов. Заметим, однако, что HMC не очень хорошо исследует положительные значения μ, что, вероятно, связано с настройками параметров шага с перешагиванием и размера шага.

c = sample(model, HMC(0.01, 10), 1000)
plot_sampler(c)
samplers 2

HMCDA

Сэмплер HMCDA является реализацией алгоритма гамильтониана Монте-Карло с двойным средним, описанным в статье The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo, авторы: Хоффман (Hoffman) и Гельман (Gelman) (2011). Интересующиеся могут найти статью на сайте arXiv.

c = sample(model, HMCDA(200, 0.65, 0.3), 1000)
plot_sampler(c)
samplers 3

MH

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

На приведенном ниже графике используется только 1000 итераций Метрополиса-Гастингса.

c = sample(model, MH(), 1000)
plot_sampler(c)
samplers 4

Как видите, сэмплер MH не очень часто перемещает оценки параметров.

NUTS

Сэмплер без разворота (No U-Turn, NUTS) является реализацией алгоритма, описанного в статье The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo, авторы: Хоффман (Hoffman) и Гельман (Gelman) (2011). Интересующиеся могут найти статью на сайте arXiv.

NUTS, как правило, очень хорошо справляется с быстрым обходом сложных апостериорных распределений.

c = sample(model, NUTS(0.65), 1000)
plot_sampler(c)
samplers 5

Единственный параметр, который необходимо задать, кроме количества итераций, — целевой коэффициент принятия. В статье Хоффмана и Гельмана отмечается, что стандартным является целевой коэффициент принятия равный 0,65.

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

c = sample(model, NUTS(0.95), 1000)
plot_sampler(c)
samplers 6

Чрезвычайно низкий коэффициент принятия будет свидетельствовать о том, что в апостериорном распределении очень мало движений:

c = sample(model, NUTS(0.2), 1000)
plot_sampler(c)
samplers 7

PG

Сэмплер Гиббса на основе частиц (PG) является реализацией алгоритма, описанного в статье Particle Markov chain Monte Carlo methods, авторы: Андрие (Andrieu), Дусе (Doucet) и Холенштайн (Holenstein) (2010). Заинтересованный читатель может найти дополнительные сведения здесь.

Два параметра — это количество частиц и количество итераций. На графике ниже показано использование 20 частиц.

c = sample(model, PG(20), 1000)
plot_sampler(c)
samplers 8

Далее построим график с 50 частицами.

c = sample(model, PG(50), 1000)
plot_sampler(c)
samplers 9