Como posso calcular uma estimativa de densidade posterior a partir de uma prévia e probabilidade?

Estou tentando entender como usar o teorema de Bayes para calcular um posterior, mas estou ficando preso à abordagem computacional, por exemplo, no seguinte caso, não está claro para mim como tirar o produto do anterior e da probabilidade e calcular o posterior:

Neste exemplo, estou interessado em calcular a probabilidade posterior de e uso um normal normal anterior em , mas quero saber como calcular o posterior a partir de um anterior em que é representado por uma cadeia MCMC, então usarei 1000 amostras como meu ponto de partida.μ p ( μ ) ∼ N ( μ = 0 , σ = 1 ) $\mu$$\mu$ $p(\mu)\sim N(\mu = 0, \sigma = 1)$$\mu$

• amostra 1000 do anterior.

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

• faça algumas observações:

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

• e calcule a probabilidade, por exemplo, :$p(y | \mu, \sigma)$

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

o que eu não entendo direito é:

1. quando / como multiplicar o anterior pela probabilidade?
2. quando / como normalizar a densidade posterior?

observe: Estou interessado na solução computacional geral que pode ser problemas generalizáveis ​​sem solução analítica

Você tem várias coisas misturadas. A teoria fala sobre multiplicar a distribuição anterior e a probabilidade, não amostras da distribuição anterior. Também não está claro do que você tem o prior, isso é um prior na média de alguma coisa? ou alguma outra coisa?

Então você tem as coisas invertidas na probabilidade, suas observações devem ser x com empates anteriores ou constantes fixas conhecidas como média e desvio padrão. E mesmo assim seria realmente o produto de 4 chamadas para dnorm com cada uma de suas observações como x e a mesma média e desvio padrão.

O que realmente não está claro é o que você está tentando fazer. Qual a sua pergunta? em quais parâmetros você está interessado? Que prévia (s) você tem desses parâmetros? existem outros parâmetros? você tem valores prévios ou fixos para eles?

Tentar fazer as coisas do jeito que você está atualmente só vai confundir você até você descobrir exatamente qual é sua pergunta e trabalhar a partir daí.

Abaixo é adicionado devido após a edição da pergunta original.

Você ainda está perdendo algumas peças e provavelmente não está entendendo tudo, mas podemos começar de onde você está.

Eu acho que você está confundindo alguns conceitos. Existe a probabilidade de mostrar a relação entre os dados e os parâmetros. Você está usando o normal, que possui 2 parâmetros, a média e o desvio padrão (ou variância ou precisão). Depois, há as distribuições anteriores nos parâmetros, você especificou um normal anterior com média 0 e sd 1, mas essa média e desvio padrão são completamente diferentes da média e desvio padrão da probabilidade. Para ser completo, você precisa conhecer a probabilidade SD ou colocar um prior na probabilidade SD, para simplificar (mas menos real) presumiremos que sabemos que a probabilidade SD é (nenhuma boa razão além de funcionar e é diferente de 1)$\frac12$

Portanto, podemos começar de maneira semelhante ao que você fez e gerar a partir do anterior:

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

Agora precisamos calcular as probabilidades, com base nos desenhos anteriores da média, na probabilidade com os dados e no valor conhecido do DP. A função dnorm nos dará a probabilidade de um único ponto, mas precisamos multiplicar os valores para cada uma das observações. Aqui está uma função para fazer isso:

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

Agora podemos calcular a probabilidade de cada sorteio do anterior para a média

> tmp <- likfun(pri)

Agora, para obter a posterior, precisamos fazer um novo tipo de sorteio, uma abordagem semelhante à amostragem por rejeição é coletar amostras dos sorteios médios anteriores proporcionais à probabilidade de cada sorteio anterior (este é o mais próximo da etapa de multiplicação que você estava perguntando sobre):

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

Agora podemos ver os resultados dos sorteios posteriores:

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

e compare os resultados acima com os valores de forma fechada da teoria

> (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

Não é uma aproximação ruim, mas provavelmente funcionaria melhor usar uma ferramenta McMC integrada para desenhar a partir do posterior. A maioria dessas ferramentas mostra um ponto de cada vez, não em lotes como acima.

Mais realisticamente, não saberíamos o SD da probabilidade e precisaríamos de um prior para isso também (geralmente o prior na variação é a ou gama), mas é mais complicado calcular (o McMC é útil) ) e não há um formulário fechado para comparar.$\chi^2$

A solução geral é usar ferramentas existentes para realizar cálculos do McMC, como WinBugs ou OpenBugs (BRugs no R fornece uma interface entre R e Bugs) ou pacotes como o LearnBayes na R.