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

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

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

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

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

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

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

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

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

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

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

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

Исчисление:

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

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

  • Ошибка SandboxViolation при установке statsmodels с помощью easy_install
  • Какой модуль статистики для python поддерживает один способ ANOVA с пост-hoc-тестами (Tukey, Scheffe или другой)?
  • Как получить прогноз теста из 2D-параметров регрессии WLS в статистических моделях
  • Построение исторической коинтеграции Значения между двумя парами
  • вычислить коэффициент определения (R2) и среднеквадратичную ошибку (RMSE) для подгонки нелинейной кривой в python
  • Python statsmodels OLS: как сохранить изученную модель в файл
  • Отсутствие перехватов регрессионных моделей OLS в статических моделях Python
  • Захват высокой многоколоничности в статистических моделях
  • 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 аргумента).

    Interesting Posts
    Python - лучший язык программирования в мире.