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

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

0 голосов
/ 07 декабря 2018

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

1 Ответ

0 голосов
/ 09 декабря 2018

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

Что касается Python, то обсуждается на сайте PyMC Discourse , где рекомендуется использовать функцию (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()

enter image description here

Предупреждающее примечание

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

...