Exemplos de amostragem Gibbs em R [fechado]

Problema

Suponha $ Y \ sim \ text {N} (\ text {média} = \ mu, \ text {Var} = \ frac {1} {\ tau}) $.

Com base em uma amostra, obtenha as distribuições posteriores de $ \ mu $ e $ \ tau $ usando o amostrador de Gibbs.

Notação

$ \ mu $ = média da população

$ \ tau $ = precisão da população (1 / variância )

$ n $ = tamanho da amostra

$ \ bar {y} $ = média da amostra

$ s ^ 2 $ = variância da amostra

Amostrador Gibbs

[ Casella, G. & George, EI (1992). Explicando o Gibbs Sampler. The American Statistician, 46, 167-174. ]

Na iteração $ i $ ($ i = 1, \ dots, N $ ):

• amostra $ \ mu ^ {(i)} $ de $ f (\ mu \, | \, \ tau ^ {(i – 1)}, \ text {data} ) $ (veja abaixo)
• amostra $ \ tau ^ {(i)} $ de $ f (\ tau \, | \, \ mu ^ {(i)}, \ text {data}) $ (veja abaixo)

A teoria garante que após um número suficientemente grande de iterações, $ T $, o conjunto $ \ {( \ mu ^ {(𝑖)}, \ tau ^ {(𝑖)}): i = T + 1, \ dots, 𝑁 \} $ pode ser visto como uma amostra aleatória da distribuição posterior conjunta.

Anteriores

$ f (\ mu, \ tau) = f (\ mu) \ vezes f (\ tau) $, com

$ f (\ mu) \ propto 1 $

$ f (\ tau) \ propto \ tau ^ {- 1} $

Posterior condicional para a média, dada a precisão $$ (\ mu \, | \, \ tau, \ text {data}) \ sim \ text {N} \ Big (\ bar {y}, \ frac {1} {n \ tau} \ Big) $$

Posterior condicional para a precisão , dada a média $$ (\ tau \, | \, \ mu, \ text {data}) \ sim \ text {Gam} \ Big (\ frac {n} {2}, \ frac {2} {(n-1) s ^ 2 + n (\ mu – \ bar {y}) ^ 2} \ Grande) $$

(rápido) Implementação R

# summary statistics of sample n <- 30 ybar <- 15 s2 <- 3 # sample from the joint posterior (mu, tau | data) mu <- rep(NA, 11000) tau <- rep(NA, 11000) T <- 1000 # burnin tau[1] <- 1 # initialisation for(i in 2:11000) { mu[i] <- rnorm(n = 1, mean = ybar, sd = sqrt(1 / (n * tau[i - 1]))) tau[i] <- rgamma(n = 1, shape = n / 2, scale = 2 / ((n - 1) * s2 + n * (mu[i] - ybar)^2)) } mu <- mu[-(1:T)] # remove burnin tau <- tau[-(1:T)] # remove burnin 

$$ $$

hist(mu) hist(tau) 

Comentários

• +1 Esta é a maneira de responder a uma pergunta aparentemente específica do código: fornecer a explicação teórica (o que a torna on-topic aqui) e então (para satisfazer o OP) forneça o código também. É bom ver você de volta!
• Obrigado pela explicação. Existe uma maneira de aprender a codificação para todos os algoritmos estatísticos. Sou um estudante e quero aprender um pouco de programação para R para poder resolver meu problema.
• Você pode me ajudar a ver esta pergunta stats.stackexchange .com / questions / 498646 / … ? Obrigado.