Интервал с наивысшей плотностью (HDI) для заднего распределения Pystan

Question

Я вижу, что в Pystan функцию HDI можно использовать для обеспечения 95% вероятного интервала, окружающего апостериорное распределение. Тем не менее, они говорят, что это будет работать только для унимодальных распределений. Если моя модель может иметь мультимодальное распределение (до 4 пиков), есть ли способ найти HDI в Pystan? Спасибо!

Answer 1

Я бы не стал рассматривать это как специфическую для Stan/PyStan проблему. Интервал с наивысшей плотностью по определению является одним интервалом и поэтому не подходит для характеристики мультимодальных распределений. Роб Хиндман (Rob Hyndman), " Вычислительные и графические регионы с наивысшей плотностью", предлагает отличную работу, которая расширяет концепцию до мультимодальных дистрибутивов, и это было реализовано в R в пакете hdrcde.

Что касается Python, это обсуждается на сайте Дискурса PyMC, где рекомендуется использовать функцию (hpd_grid) что Освальдо Мартин написал для своей книги "Байесовский анализ с Python". Источник для этой функции находится в файле hpd.py, а для области 95% будет использоваться как

from hpd import hpd_grid

intervals, x, y, modes = hpd_grid(samples, alpha=0.05)

где samples являются задними образцами одного из ваших параметров, и intervals список кортежей, представляющих регионы с наибольшей плотностью

Пример с участком

Вот пример графика с использованием поддельных мультимодальных данных.

import numpy as np
from matplotlib import pyplot as plt
from hpd import hpd_grid

# include two modes
samples = np.random.normal(loc=[-4,4], size=(1000, 2)).flatten()

# compute high density regions
hpd_mu, x_mu, y_mu, modes_mu = hpd_grid(samples)

plt.figure(figsize=(8,6))

# raw data
plt.hist(samples), density=True, bins=29, alpha=0.5)

# estimated distribution
plt.plot(x_mu, y_mu)

# high density intervals
for (x0, x1) in hpd_mu:
    plt.hlines(y=0, xmin=x0, xmax=x1, linewidth=5)
    plt.axvline(x=x0, color='grey', linestyle='--', linewidth=1)
    plt.axvline(x=x1, color='grey', linestyle='--', linewidth=1)

# modes
for xm in modes_mu:
    plt.axvline(x=xm, color='r')

plt.show()

---

Предупреждение

Следует отметить, что мультимодальные апостериорные распределения по правильно смоделированным параметрам, как правило, встречаются редко, но встречаются очень часто при неконвергентной выборке MCMC, особенно при использовании нескольких цепочек (что является наилучшей практикой). Если кто-то ожидает мультимодальность априори, то обычно это приводит к некоторой форме смешанной модели, которая устранит мультимодальность. Если вы не ожидаете мультимодальности, но постеры демонстрируют это в любом случае, то это красный флаг, который скептически относится к результатам.