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

С учетом заднего p (Θ | D) по некоторым параметрам Θ можно определить следующее:

Самая высокая плотность залегания:

Самая высокая по задняя плотность – это набор наиболее вероятных значений Θ, которые в совокупности составляют 100 (1-α)% от задней массы.

Другими словами, для данного α мы ищем a * p **, который удовлетворяет:

введите описание изображения здесь

а затем получить самый высокий уровень плотности задних частот в качестве набора:

введите описание изображения здесь

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

Используя те же обозначения, что и выше, заслуживающий доверия регион (или интервал) определяется как:

введите описание изображения здесь

В зависимости от распределения может быть много таких интервалов. Центральный надежный интервал определяется как надежный интервал, где на каждом хвосте имеется (1-α) / 2 масса.

Исчисление:

  • Для общих распределений, данных выборок из дистрибутива, есть ли какие-либо встроенные модули для получения двух величин выше в Python или PyMC ?

  • Для общих параметрических распределений (например, Beta, Gaussian и т. Д.) Существуют ли какие-либо встроенные модули или библиотеки для вычисления этого с помощью SciPy или statsmodels ?

6 Solutions collect form web for “Наивысшая область засушливой плотности и Центральный заслуживающий доверия регион”

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

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

По моему мнению, «центральный заслуживающий доверия регион» ничем не отличается от того, как рассчитываются доверительные интервалы; все, что вам нужно, это инверсия функции 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, как p* = 0.0450 ниже; красная область представляет 1 - alpha распределения, а горизонтальная пунктирная линия – p* .

введите описание изображения здесь

PyMC имеет встроенную функцию для вычисления hpd. В версии 2.3 он находится в утилях. См. Источник здесь . В качестве примера линейной модели и 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% -ный центральный вероятный интервал.

Другой вариант (адаптированный от 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 от scipy.

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

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

является

  [-1.95996398 1.95996398] 

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

Я наткнулся на это сообщение, пытаясь найти способ оценить ИРЧП из образца 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) 

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

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

  • Как добавить функции регрессии в python или создать новую функцию регрессии из заданных коэффициентов?
  • Какой модуль статистики для python поддерживает один способ ANOVA с пост-hoc-тестами (Tukey, Scheffe или другой)?
  • Сравнение результатов от статмоделей ARIMA с исходными данными
  • Как перебирать столбцы базы данных pandas для запуска регрессии
  • statsmodels linear regression - patsy formula для включения всех предикторов в модель
  • Каковы подводные камни использования Dill для сериализации моделей scikit-learn / statsmodels?
  • Вычислить логическую регрессию в python
  • Python - Прокручивающееся окно Оценка регрессии OLS
  • Python - лучший язык программирования в мире.