Наивысшая область засушливой плотности и Центральный заслуживающий доверия регион

Ответ 1

Чтобы вычислить HPD, вы можете использовать pymc3. Вот пример

import pymc3
from scipy.stats import norm
a = norm.rvs(size=10000)
pymc3.stats.hpd(a)

Ответ 2

По моему мнению, "центральный заслуживающий доверия регион" не отличается от того, как рассчитываются доверительные интервалы; все, что вам нужно, это инверсия функции cdf в alpha/2 и 1-alpha/2; в scipy это называется ppf (функция процентных точек); так как для гауссовского заднего распределения:

>>> from scipy.stats import norm
>>> alpha = .05
>>> l, u = norm.ppf(alpha / 2), norm.ppf(1 - alpha / 2)

чтобы убедиться, что [l, u] покрывает (1-alpha) задней плотности:

>>> norm.cdf(u) - norm.cdf(l)
0.94999999999999996

аналогично для Beta posterior, скажем a=1 и b=3:

>>> from scipy.stats import beta
>>> l, u = beta.ppf(alpha / 2, a=1, b=3), beta.ppf(1 - alpha / 2, a=1, b=3)

и снова:

>>> beta.cdf(u, a=1, b=3) - beta.cdf(l, a=1, b=3)
0.94999999999999996

здесь вы можете увидеть параметрические распределения, включенные в scipy; и я думаю, что у всех из них есть функция ppf;

Что касается самой высокой области задней плотности, это более сложно, так как функция pdf не обязательно обратима; и вообще такая область может быть даже не связана; например, в случае бета с a = b = .5 (как можно видеть здесь);

Но, в случае распределения Гаусса, нетрудно видеть, что "область с наивысшей плотностью залегания" совпадает с "центральным заслуживающим доверия регионом"; и я думаю, что это имеет место для всех симметричных унимодальных распределений (т.е. если функция PDF симметрична относительно способа распределения)

Возможным численным подходом для общего случая будет бинарный поиск по значению p* с использованием численного интегрирования pdf; используя тот факт, что интеграл является монотонной функцией от p*;

Вот пример смеси Гаусса:

[1] Прежде всего вам нужна аналитическая функция pdf; для смеси Гаусса легко:

def mix_norm_pdf(x, loc, scale, weight):
    from scipy.stats import norm
    return np.dot(weight, norm.pdf(x, loc, scale))

так, например, для значений местоположения, масштаба и веса, как в

loc    = np.array([-1, 3])   # mean values
scale  = np.array([.5, .8])  # standard deviations
weight = np.array([.4, .6])  # mixture probabilities

вы получите два хороших распределения Гаусса, держась за руки:

[2] теперь вам нужна функция ошибки, которая задает тестовое значение для p*, интегрирует функцию PDF выше p* и возвращает квадратную ошибку от желаемого значения 1 - alpha:

def errfn( p, alpha, *args):
    from scipy import integrate
    def fn( x ):
        pdf = mix_norm_pdf(x, *args)
        return pdf if pdf > p else 0

    # ideally integration limits should not
    # be hard coded but inferred
    lb, ub = -3, 6 
    prob = integrate.quad(fn, lb, ub)[0]
    return (prob + alpha - 1.0)**2

[3] теперь при заданном значении alpha мы можем минимизировать функцию ошибки, чтобы получить p*:

alpha = .05

from scipy.optimize import fmin
p = fmin(errfn, x0=0, args=(alpha, loc, scale, weight))[0]

что приводит к p* = 0.0450 и HPD, как показано ниже; красная область представляет 1 - alpha распределения, а горизонтальная пунктирная линия p*.

Ответ 3

PyMC имеет встроенную функцию для вычисления hpd. В версии 3.2 она используется. См. Источник здесь. В качестве примера линейной модели и HPD

import pymc as pc  
import numpy as np
import matplotlib.pyplot as plt 
## data
np.random.seed(1)
x = np.array(range(0,50))
y = np.random.uniform(low=0.0, high=40.0, size=50)
y = 2*x+y
## plt.scatter(x,y)

## priors
emm = pc.Uniform('m', -100.0, 100.0, value=0)
cee = pc.Uniform('c', -100.0, 100.0, value=0) 

#linear-model
@pc.deterministic(plot=False)
def lin_mod(x=x, cee=cee, emm=emm):
    return emm*x + cee 

#likelihood
llhy = pc.Normal('y', mu=lin_mod, tau=1.0/(10.0**2), value=y, observed=True)

linearModel = pc.Model( [llhy, lin_mod, emm, cee] )
MCMClinear = pc.MCMC( linearModel)
MCMClinear.sample(10000,burn=5000,thin=5)
linear_output=MCMClinear.stats()

## pc.Matplot.plot(MCMClinear)
## print HPD using the trace of each parameter 
print(pc.utils.hpd(MCMClinear.trace('m')[:] , 1.- 0.95))
print(pc.utils.hpd(MCMClinear.trace('c')[:] , 1.- 0.95))

Вы также можете рассмотреть возможность расчета квантилей

print(linear_output['m']['quantiles'])
print(linear_output['c']['quantiles'])

где я думаю, что если вы просто возьмете 2,5% до 97,5%, вы получите свой 95% -ный центральный надежный интервал.

Ответ 4

Другой вариант (адаптированный от R к Python) и взятый из книги Джона К. Крушке "Анализ байесовских данных" заключается в следующем:

from scipy.optimize import fmin
from scipy.stats import *

def HDIofICDF(dist_name, credMass=0.95, **args):
    # freeze distribution with given arguments
    distri = dist_name(**args)
    # initial guess for HDIlowTailPr
    incredMass =  1.0 - credMass

    def intervalWidth(lowTailPr):
        return distri.ppf(credMass + lowTailPr) - distri.ppf(lowTailPr)

    # find lowTailPr that minimizes intervalWidth
    HDIlowTailPr = fmin(intervalWidth, incredMass, ftol=1e-8, disp=False)[0]
    # return interval as array([low, high])
    return distri.ppf([HDIlowTailPr, credMass + HDIlowTailPr])

Идея состоит в том, чтобы создать функцию intervalWidth, которая возвращает ширину интервала, который начинается с lowTailPr и имеет массу credMass. Минимум функции intervalWidth основан на использовании Fmin Minimizer из SciPy.

Например, результат:

print HDIofICDF(norm, credMass=0.95, loc=0, scale=1)

является

    [-1.95996398  1.95996398]

Имя параметров распространения, передаваемых в HDIofICDF, должно быть точно таким же, как и в scipy.

Ответ 5

Я наткнулся на это сообщение, пытаясь найти способ оценить ИРЧП из образца MCMC, но ни один из ответов не работал у меня. Как aloctavodia, я адаптировал пример R из книги Doing Bayesian Data Analysis для Python. Мне нужно было вычислить 95% ИРЧП из образца MCMC. Здесь мое решение:

import numpy as np
def HDI_from_MCMC(posterior_samples, credible_mass):
    # Computes highest density interval from a sample of representative values,
    # estimated as the shortest credible interval
    # Takes Arguments posterior_samples (samples from posterior) and credible mass (normally .95)
    sorted_points = sorted(posterior_samples)
    ciIdxInc = np.ceil(credible_mass * len(sorted_points)).astype('int')
    nCIs = len(sorted_points) - ciIdxInc
    ciWidth = [0]*nCIs
    for i in range(0, nCIs):
    ciWidth[i] = sorted_points[i + ciIdxInc] - sorted_points[i]
    HDImin = sorted_points[ciWidth.index(min(ciWidth))]
    HDImax = sorted_points[ciWidth.index(min(ciWidth))+ciIdxInc]
    return(HDImin, HDImax)

Метод выше дает мне логические ответы на основе данных, которые у меня есть!

Ответ 6

Вы можете получить центральный надежный интервал двумя способами: графически, когда вы вызываете summary_plot для переменных в вашей модели, по умолчанию установлен флаг bpd, который установлен на True. Изменение этого параметра на False будет отображать центральные интервалы. Второе место вы можете получить, когда вы вызываете метод summary на вашей модели или node; он даст вам задние квантиля, а внешние будут по умолчанию 95% -ным центральным интервалом (который вы можете изменить с помощью аргумента alpha).