Kommentare
- Ich habe ein Buch mit dem Titel geschrieben, in dem Monte-Carlo-Methoden mit R vorgestellt werden, die Sie überprüfen können für solche Beispiele.
Antwort
Problem
Angenommen, $ Y \ sim \ text {N} (\ text {mean} = \ mu, \ text {Var} = \ frac {1} {\ tau}) $.
Ermitteln Sie anhand einer Stichprobe die posterioren Verteilungen von $ \ mu $ und $ \ tau $ mit dem Gibbs-Sampler.
Notation
$ \ mu $ = Populationsmittel
$ \ tau $ = Populationsgenauigkeit (1 / Varianz )
$ n $ = Stichprobengröße
$ \ bar {y} $ = Stichprobenmittelwert
$ s ^ 2 $ = Stichprobenvarianz
Gibbs-Sampler
[ Casella, G. & George, EI (1992). Erklären des Gibbs-Samplers. The American Statistician, 46, 167–174. ]
Bei Iteration $ i $ ($ i = 1, \ dots, N $ ):
- Beispiel $ \ mu ^ {(i)} $ aus $ f (\ mu \, | \, \ tau ^ {(i – 1)}, \ text {data} ) $ (siehe unten)
- Beispiel $ \ tau ^ {(i)} $ aus $ f (\ tau \, | \, \ mu ^ {(i)}, \ text {data}) $ (siehe unten)
Die Theorie stellt sicher, dass nach einer ausreichend großen Anzahl von Iterationen $ T $ die Menge $ \ {( \ mu ^ {(𝑖)}, \ tau ^ {(𝑖)}): i = T + 1, \ dots, 𝑁 \} $ kann als Zufallsstichprobe aus der gemeinsamen posterioren Verteilung angesehen werden.
Priors
$ f (\ mu, \ tau) = f (\ mu) \ mal f (\ tau) $, mit
$ f (\ mu) \ propto 1 $
$ f (\ tau) \ propto \ tau ^ {- 1} $
Bedingter posteriorer Wert für den Mittelwert angesichts der Genauigkeit $$ (\ mu \, | \, \ tau, \ text {data}) \ sim \ text {N} \ Big (\ bar {y}, \ frac {1} {n \ tau} \ Big) $$
Bedingter posteriorer Wert für die Präzision unter Berücksichtigung des Mittelwerts $$ (\ tau \, | \, \ mu, \ text {data}) \ sim \ text {Gam} \ Big (\ frac {n}) {2}, \ frac {2} {(n-1) s ^ 2 + n (\ mu – \ bar {y}) ^ 2} \ Big) $$
(schnelle) R-Implementierung
# 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) 
Kommentare
- +1 Dies ist der Weg, um eine scheinbar codespezifische Frage zu beantworten: Geben Sie die theoretische Erklärung (die es hier zum Thema macht) und Geben Sie dann (um das OP zu erfüllen) auch den Code an. Schön, Sie wiederzusehen!
- Vielen Dank für die Erklärung. Gibt es eine Möglichkeit, die Codierung für alle statistischen Algorithmen zu lernen? Ich bin Student und möchte etwas Codierung für R lernen, damit ich mein Problem lösen kann.
- Können Sie mir helfen, diese Frage zu sehen? stats.stackexchange .com / question / 498646 / … ? Vielen Dank.