¿Cómo puedo calcular una estimación de densidad posterior a partir de una probabilidad previa y?

Question

Estoy tratando de entender cómo utilizar el teorema de Bayes para calcular la parte posterior, pero estoy pegado con el método de cálculo, por ejemplo, en el siguiente caso no es claro para mí cómo tomar el producto de la previa y la probabilidad y, a continuación, calcular la parte posterior:

Para este ejemplo, estoy interesado en calcular la probabilidad posterior de $\mu$, y el uso de un estándar normal antes de $\mu$ $p(\mu)\sim N(\mu = 0, \sigma = 1)$, pero quiero saber cómo calcular la parte posterior de un anterior en $\mu$ que es representado por un MCMC de la cadena, así que voy a usar 1000 muestras como mi punto de partida.

• muestra de 1000 de los anteriores.

  set.seed(0)
  prior.mu      <- 0
  prior.sigma   <- 1
  prior.samples <- sort(rnorm(1000, prior.mu, prior.sigma))
  

• hacer algunas observaciones:

  observations <- c(0.4, 0.5, 0.8, 0.1)
  

• y calcular la probabilidad, por ejemplo,$p(y | \mu, \sigma)$:

  likelihood <- prod(dnorm(observations, mean(prior.samplse), sd(prior.samples)))
  

lo que no acabo de entender es:

1. cuando / cómo multiplicar la anterior por la probabilidad?
2. cuando / cómo normalizar la parte posterior de la densidad?

por favor nota: estoy interesado en la general computacional de la solución, que podría ser generalizable problemas sin solución analítica

Answer 1

Tiene varias cosas mezcladas. La teoría habla de la multiplicación de la previa de la distribución y de la probabilidad, no de muestras de la distribución previa. También no está claro lo que usted tiene antes de, esta es una previa de la media de algo? o algo más?

Entonces usted tiene cosas invertido en la probabilidad, sus observaciones debe ser x con cualquiera de los anteriores sorteos o fijo conocido constantes como la media y la desviación estándar. E incluso entonces realmente sería el producto de 4 llamadas a dnorm con cada una de sus observaciones como x y la misma media y desviación estándar.

Lo que realmente no es muy claro es lo que usted está tratando de hacer. ¿Cuál es su pregunta? qué parámetros son los que te interesan? lo que antes de la(s) tiene en esos parámetros? hay otros parámetros? ¿tienes dudas o valores fijos para los?

Tratando de ir sobre las cosas de la manera que en la actualidad son sólo va a confundir más hasta que el trabajo fuera exactamente lo que su pregunta es y trabajar desde allí.

A continuación se añade vence después de la edición de la pregunta original.

Aún le faltan algunas piezas, y probablemente no entienda todo, pero podemos empezar desde donde estás.

Creo que se están confundiendo algunos conceptos. Existe la probabilidad de que se muestra la relación entre los datos y los parámetros, que están utilizando la normal que tiene 2 parámetros, la media y la desviación estándar (o la varianza, o de precisión). Luego está el antes de las distribuciones de los parámetros, se han especificado previa normal con media 0 y sd 1, pero que la media y la desviación estándar son completamente diferentes a los de la media y la desviación estándar de la probabilidad. Para ser completa debe saber la probabilidad SD o colocar un previo acerca de la posibilidad de SD, por simplicidad (pero menos real) voy a suponer que sabemos que la probabilidad de SD es $\frac12$ (ninguna buena razón por la que no funciona y que es diferente de 1).

Así que podemos empezar similar a lo que hizo y generar a partir de la anterior:

> obs <- c(0.4, 0.5, 0.8, 0.1)
> pri <- rnorm(10000, 0, 1)

Ahora tenemos que calcular las probabilidades, esto se basa en los anteriores sorteos de la media, la probabilidad con los datos, y el valor conocido de la SD. El dnorm función nos dará la probabilidad de que un solo punto, pero tenemos que multiplicar los valores para cada una de las observaciones, aquí es una función para hacer que:

> likfun <- function(theta) {
+ sapply( theta, function(t) prod( dnorm(obs, t, 0.5) ) )
+ }

Ahora podemos calcular la probabilidad de cada sorteo de la previa de la media

> tmp <- likfun(pri)

Ahora para obtener la parte posterior tenemos que hacer un nuevo tipo de sorteo, un enfoque que es similar a la de rechazo de muestreo es la muestra de la antes de la media dibuja proporcional a la probabilidad de cada uno antes del sorteo (esta es la más cercana a la multiplicación de paso le preguntaba acerca de):

> post <- sample( pri, 100000, replace=TRUE, prob=tmp )

Ahora podemos ver los resultados de la parte posterior de los sorteos:

> mean(post)
[1] 0.4205842
> sd(post)
[1] 0.2421079
> 
> hist(post)
> abline(v=mean(post), col='green')

y comparar los resultados anteriores a la forma cerrada de valores a partir de la teoría

> (1/1^2*mean(pri) + length(obs)/0.5^2 * mean(obs))/( 1/1^2 + length(obs)/0.5^2 )
[1] 0.4233263
> sqrt(1/(1+4*4))
[1] 0.2425356

No es una mala aproximación, pero es probable que funcione mejor para usar un built-in McMC herramienta para dibujar de la parte posterior. La mayoría de estas herramientas se muestra un punto en un momento en lotes de arriba.

De manera más realista que no sé la SD de la probabilidad y sería necesario un previo para que así (a menudo antes de la varianza es una $\chi^2$ o gamma), pero luego es más complicado de calcular (McMC viene muy bien) y no hay forma cerrada para comparar.

La solución general es utilizar las herramientas existentes para hacer la McMC cálculos tales como WinBugs o OpenBugs (BRugs en R proporciona una interfaz entre la investigación y Errores) o paquetes de LearnBayes en R.