importar - statsmodels python
Región de mayor densidad posterior y región central creíble (6)
A mi entender, "región central creíble" no es diferente de cómo se calculan los intervalos de confianza; todo lo que necesita es el inverso de la función cdf en alpha/2 y 1-alpha/2 ; en scipy esto se llama ppf (función de punto porcentual); así como para la distribución posterior gaussiana:
>>> from scipy.stats import norm
>>> alpha = .05
>>> l, u = norm.ppf(alpha / 2), norm.ppf(1 - alpha / 2)
para verificar que [l, u] cubre (1-alpha) de densidad posterior:
>>> norm.cdf(u) - norm.cdf(l)
0.94999999999999996
de manera similar para Beta posterior con, por ejemplo, a=1 y 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)
y otra vez:
>>> beta.cdf(u, a=1, b=3) - beta.cdf(l, a=1, b=3)
0.94999999999999996
here puede ver las distribuciones paramétricas que se incluyen en scipy; y supongo que todos ellos tienen función ppf ;
En cuanto a la región de densidad posterior más alta, es más difícil, ya que la función pdf no es necesariamente invertible; y en general tal región ni siquiera puede estar conectada; por ejemplo, en el caso de Beta con a = b = .5 (como se puede ver here );
Pero, en el caso de la distribución gaussiana, es fácil ver que la "región de densidad posterior más alta" coincide con la "región central creíble"; y creo que ese es el caso de todas las distribuciones uni-modales simétricas (es decir, si la función pdf es simétrica alrededor del modo de distribución)
Un posible enfoque numérico para el caso general sería la búsqueda binaria sobre el valor de p* utilizando la integración numérica de pdf ; utilizando el hecho de que la integral es una función monótona de p* ;
Aquí hay un ejemplo de mezcla gaussiana:
[1] Lo primero que necesitas es una función de pdf analítica; Para la mezcla gaussiana que es fácil:
def mix_norm_pdf(x, loc, scale, weight):
from scipy.stats import norm
return np.dot(weight, norm.pdf(x, loc, scale))
así, por ejemplo, para los valores de ubicación, escala y peso como en
loc = np.array([-1, 3]) # mean values
scale = np.array([.5, .8]) # standard deviations
weight = np.array([.4, .6]) # mixture probabilities
obtendrás dos agradables distribuciones gaussianas cogidas de la mano:
[2] ahora, necesita una función de error que dado un valor de prueba para p* integra la función pdf por encima de p* y devuelve un error cuadrado del valor deseado 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] ahora, para un valor dado de alpha podemos minimizar la función de error para obtener p* :
alpha = .05
from scipy.optimize import fmin
p = fmin(errfn, x0=0, args=(alpha, loc, scale, weight))[0]
que da como resultado p* = 0.0450 y HPD como se muestra a continuación; el área roja representa 1 - alpha de la distribución, y la línea discontinua horizontal es p* .
Dada una p posterior (Θ | D) sobre algunos parámetros Θ, se puede define lo siguiente:
Región de mayor densidad posterior:
La región de densidad posterior más alta es el conjunto de valores más probables de que, en total, constituyen el 100 (1-α)% de la masa posterior.
En otras palabras, para un α dado, buscamos un * p ** que satisfaga:
y luego obtener la región de mayor densidad posterior como el conjunto:
Región central creíble:
Usando la misma notación que arriba, una Región creíble (o intervalo) se define como:
Dependiendo de la distribución, podría haber muchos tales intervalos. El intervalo creíble central se define como un intervalo creíble donde hay una masa (1-α) / 2 en cada cola .
Cálculo:
Para las distribuciones generales, dadas las muestras de la distribución, ¿hay algunas incorporaciones para obtener las dos cantidades anteriores en Python o PyMC ?
Para distribuciones paramétricas comunes (por ejemplo, Beta, Gaussian, etc.), ¿hay algunas incorporadas o bibliotecas para calcular esto usando SciPy o statsmodels ?
Me topé con esta publicación tratando de encontrar una manera de estimar un IDH a partir de una muestra de MCMC pero ninguna de las respuestas funcionó para mí. Como aloctavodia, adapté un ejemplo de R del libro Doing Bayesian Data Analysis a Python. Necesitaba calcular un HDI del 95% a partir de una muestra de MCMC. Aquí está mi solución:
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)
¡El método anterior me está dando respuestas lógicas basadas en los datos que tengo!
Otra opción (adaptada de R a Python) y tomada del libro Análisis de datos bayesianos de John K. Kruschke) es la siguiente:
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])
La idea es crear una función intervalWidth que devuelva el ancho del intervalo que comienza en lowTailPr y tiene masa de credMass . El mínimo de la función intervalWidth se basa utilizando el minimizador fmin desde scipy.
Por ejemplo el resultado de:
print HDIofICDF(norm, credMass=0.95, loc=0, scale=1)
es
[-1.95996398 1.95996398]
El nombre de los parámetros de distribución pasados a HDIofICDF, debe ser exactamente el mismo utilizado en scipy.
Para calcular HPD puede aprovechar pymc3, Aquí hay un ejemplo
import pymc3
from scipy.stats import norm
a = norm.rvs(size=10000)
pymc3.stats.hpd(a)
Puede obtener el intervalo central creíble de dos maneras: gráficamente, cuando llama a summary_plot en las variables de su modelo, hay un indicador de bpd que se establece en True de manera predeterminada. Cambiar esto a False dibujará los intervalos centrales. El segundo lugar donde puede obtenerlo es cuando llama al método de summary en su modelo o nodo; le dará cuantiles posteriores, y los externos serán un intervalo central del 95% por defecto (que puede cambiar con el argumento alpha ).
PyMC tiene una función incorporada para calcular el hpd. En v2.3 está en utils. Vea la fuente here . Como ejemplo de un modelo lineal y su 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))
También puedes considerar calcular los cuantiles
print(linear_output[''m''][''quantiles''])
print(linear_output[''c''][''quantiles''])
donde creo que si solo tomas los valores de 2.5% a 97.5% obtendrás tu intervalo creíble central de 95%.