Inferir la distribución previa

Question

Supongamos que tomamos una muestra ( $X_1, X_2, ... X_n$ ) de una distribución en la que suponemos que $X_i $ ~ $ Bin(n_i, p_i)$ y $n_i$ es conocido por cada $i$ . También suponemos que $p_i$ son independientes e idénticamente distribuidos, $p_i$ ~ $D$ , donde $D$ es una distribución desconocida. $n_i$ no puede suponerse que sea grande.

Mi objetivo es obtener una estimación bayesiana (o una distribución de probabilidad) para $p_i$ . Pero para ello es necesario idear una distribución para $D$ .

Una opción es hacer una distribución empírica que utilice estimaciones frecuentistas para cada $p_i$ (es decir $p_i = X_i/n_i$ ). Esta es una idea bastante intuitiva y potencialmente razonable. Por desgracia, la presencia de pequeños $n_i$ haría que las colas fueran más pesadas de lo que deberían (muchos valores extremos cercanos a 0 o 1).

Estoy buscando otra opción que no tenga los problemas de la solución mencionada.

Una posibilidad que tengo en mente es utilizar el siguiente algoritmo:

1. Generar la distribución previa como se ha explicado anteriormente.
2. Obtenga una estimación de MAP o EAP para cada $p_i$ .
3. Generar una nueva prioridad empírica a partir de las probabilidades obtenidas en 2.
4. Vuelve a 2 (continúa durante un número determinado de pasos, o posiblemente hasta la convergencia )

¿Es este método similar a cualquier otro que exista? ¿Es razonable?

Answer 1

¡Espero que te guste Python! Voy a recitar mi comentario aquí:

Esto parece un modelo jerárquico. Si quisiera recrear el conjunto de datos, esto es lo que haría: Que $D$ ser un $Beta(\alpha, \beta)$ distribución (razonable ya que se trata de probabilidades). No sabemos $\alpha, \beta$ les asignamos priores, digamos exponenciales para ambos con algún $\lambda$ hiperparámetro. A continuación, dibujamos el $p_i$ para cada $i$ y muestra $X_i$ de los binomios.

Así es como recrearía el conjunto de datos. Para hacer inferencia, vamos hacia atrás. Aquí está el modelo en PyMC:

import pymc as pm

#fake data
X = np.array([3,2,2,5,7,10,11])
n = np.array([5, 4, 4, 6, 10, 19, 12])

#here I make sure I fulfill fake-data constraints
assert X.shape == n.shape
assert (X <= n).all()

alpha = pm.Exponential("alpha", 1)
beta = pm.Exponential("beta", 1)

p = pm.Beta( "p", alpha, beta, size=X.shape[0])

obs = pm.Binomial("obs", n, p, value=X, observed=True)

mcmc = pm.MCMC([obs,p,beta,alpha])
mcmc.sample(10000, 5000)

Y alguna salida:

Con muestras de las posteriors de $\alpha$ y $\beta$ podemos reconstruir las posibles distribuciones de $D$ la distribución desconocida:

Editar : Disculpas, el eje x debería estar entre 0-1, no 0-500, es algo de Python que olvidé cambiar.

Answer 2

El algoritmo que describe está tratando $D$ como una variable en el problema, pero utilizando un método distinto a la inferencia bayesiana para tratar $D$ . Una solución bayesiana requiere manejar $D$ de forma bayesiana, es decir, dándole una prioridad e integrándola para obtener el marginal de $p_i$ . Por ejemplo, se puede utilizar un proceso Dirichlet como la prioridad para $D$ y hacer inferencia para $p_i$ mediante el muestreo de Gibbs.