Визуализация сэмплера
Код
Для каждого сэмплера мы будем использовать один и тот же код для построения графиков его путей. Приведенный ниже блок загружает соответствующие библиотеки и определяет функцию для построения графика траектории сэмплера по апостериорному распределению.
Используемое здесь определение модели 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)
HMC
Выборка гамильтониана Монте-Карло (HMC) — это типичный используемый сэмплер, поскольку он достаточно хорошо сходится. Однако часто бывает сложно задать правильные параметры для этого сэмплера, и сэмплер NUTS часто проще запустить, если вы не хотите тратить много времени на определение размера шага и количества шагов. Заметим, однако, что HMC не очень хорошо исследует положительные значения μ, что, вероятно, связано с настройками параметров шага с перешагиванием и размера шага.
c = sample(model, HMC(0.01, 10), 1000)
plot_sampler(c)
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)
MH
Выборка Метрополиса-Гастингса (MH) представляет собой один из самых ранних методов Монте-Карло с цепью Маркова. Выборка MH выполняет не так много проходов, в отличие от многих других сэмплеров, реализованных в Turing. Как правило, для сходимости к подходящей оценке параметров требуется гораздо более длинная цепочка.
На приведенном ниже графике используется только 1000 итераций Метрополиса-Гастингса.
c = sample(model, MH(), 1000)
plot_sampler(c)
Как видите, сэмплер 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)
Единственный параметр, который необходимо задать, кроме количества итераций, — целевой коэффициент принятия. В статье Хоффмана и Гельмана отмечается, что стандартным является целевой коэффициент принятия равный 0,65.
Вот график, показывающий очень высокий коэффициент принятия. Обратите внимание, что он, похоже, зависит от режима и не особенно подходит для исследования апостериорного распределения по сравнению с вариантом целевого коэффициента принятия, равного 0,65.
c = sample(model, NUTS(0.95), 1000)
plot_sampler(c)
Чрезвычайно низкий коэффициент принятия будет свидетельствовать о том, что в апостериорном распределении очень мало движений:
c = sample(model, NUTS(0.2), 1000)
plot_sampler(c)
PG
Сэмплер Гиббса на основе частиц (PG) является реализацией алгоритма, описанного в статье Particle Markov chain Monte Carlo methods, авторы: Андрие (Andrieu), Дусе (Doucet) и Холенштайн (Holenstein) (2010). Заинтересованный читатель может найти дополнительные сведения здесь.
Два параметра — это количество частиц и количество итераций. На графике ниже показано использование 20 частиц.
c = sample(model, PG(20), 1000)
plot_sampler(c)
Далее построим график с 50 частицами.
c = sample(model, PG(50), 1000)
plot_sampler(c)