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

Estoy tratando de entender cómo usar el teorema de Bayes para calcular un posterior, pero me estoy atascando con el enfoque computacional, por ejemplo, en el siguiente caso no me queda claro cómo tomar el producto de lo anterior y la probabilidad y luego calcular el posterior:

Para este ejemplo, estoy interesado en calcular la probabilidad posterior de y utilizo un estándar normal previo en \ mu p (\ mu) \ sim N (\ mu = 0, \ sigma = 1) , pero quiero saber cómo calcular la parte posterior de una anterior en \ mu que está representada por una cadena MCMC, por lo que utilizaré 1000 muestras como punto de partida.$\mu$$\mu$ $p(\mu)\sim N(\mu = 0, \sigma = 1)$$\mu$

• muestra 1000 de la anterior.

  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 calcule la probabilidad, por ejemplo, $p(y | \mu, \sigma)$ :

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

lo que no entiendo es:

1. ¿Cuándo / cómo multiplicar lo anterior por la probabilidad?
2. ¿Cuándo / cómo normalizar la densidad posterior?

tenga en cuenta: estoy interesado en la solución computacional general que podría ser problemas generalizables sin solución analítica

Tienes varias cosas mezcladas. La teoría habla de multiplicar la distribución previa y la probabilidad, no muestras de la distribución previa. Además, no está claro de qué tienes el prior, ¿es un prior en el sentido de algo? ¿o algo mas?

Entonces tiene las cosas invertidas en la probabilidad, sus observaciones deben ser x con dibujos anteriores o constantes fijas conocidas como la media y la desviación estándar. Y aun así, 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 está claro es lo que está tratando de hacer. ¿Cuál es tu pregunta? ¿Qué parámetros le interesan? ¿Qué prior (es) tienes sobre esos parámetros? ¿hay otros parámetros? ¿tiene valores previos o fijos para esos?

Tratar de hacer las cosas de la forma en que se encuentra actualmente solo lo confundirá más hasta que resuelva exactamente cuál es su pregunta y trabaje desde allí.

A continuación se agrega después de la edición de la pregunta original.

Todavía te faltan algunas piezas y probablemente no entiendas todo, pero podemos comenzar desde donde estás.

Creo que estás confundiendo algunos conceptos. Existe la probabilidad de que muestre la relación entre los datos y los parámetros, está utilizando la normal que tiene 2 parámetros, la media y la desviación estándar (o varianza o precisión). Luego están las distribuciones previas en los parámetros, usted ha especificado un previo normal con media 0 y sd 1, pero esa media y la desviación estándar son completamente diferentes de la media y la desviación estándar de la probabilidad. Para completar, debe conocer la probabilidad de SD o colocar un previo en la probabilidad de SD, por simplicidad (pero menos real) Asumiré que sabemos que la probabilidad de SD es (no hay una buena razón más que funciona y es diferente de 1)$\frac12$

Entonces podemos comenzar de manera similar a lo que hizo y generar a partir de lo 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 dibujos anteriores de la media, la probabilidad con los datos y el valor conocido de la SD. La función dnorm nos dará la probabilidad de un solo punto, pero necesitamos multiplicar los valores para cada una de las observaciones, aquí hay una función para hacer eso:

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

Ahora podemos calcular la probabilidad para cada sorteo del anterior para la media

> tmp <- likfun(pri)

Ahora para obtener el posterior necesitamos hacer un nuevo tipo de sorteo, un enfoque que es similar al muestreo de rechazo es tomar muestras de los sorteos medios anteriores proporcionales a la probabilidad de cada sorteo anterior (este es el paso más cercano al de multiplicación preguntar por):

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

Ahora podemos ver los resultados de los sorteos posteriores:

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

y compare los resultados anteriores con los valores de forma cerrada 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 probablemente funcionaría mejor usar una herramienta McMC incorporada para dibujar desde la parte posterior. La mayoría de estas herramientas muestra un punto a la vez, no en lotes como el anterior.

De manera más realista, no sabríamos la SD de la probabilidad y también necesitaríamos un previo para eso (a menudo, el previo en la varianza es un o gamma), pero luego es más complicado de calcular (McMC es útil ) y no hay una forma cerrada para comparar.$\chi^2$

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