To pytanie jest techniczną kontynuacją tego pytania .
Mam problem ze zrozumieniem i replikacją modelu przedstawionego w Raftery (1988): Wnioskowanie o dwumianowym parametrze : hierarchiczne podejście Bayesa w WinBUGS / OpenBUGS / JAGS. Nie chodzi tylko o kod, więc powinien być tutaj na temat.
tło
Niech będzie zbiorem liczb sukcesów z rozkładu dwumianowego o nieznanych i . Ponadto zakładam, że ma rozkład Poissona z parametrem (jak omówiono w artykule). Następnie każdy ma rozkład Poissona ze średnią . Chcę określić priorytety w kategoriach i .N θ N μ x i λ = μ θ λ θ
Zakładając, że nie mam żadnej dobrej wcześniejszej wiedzy na temat lub , chcę przypisać nieinformacyjne priory zarówno i . Powiedzmy, że moimi priorytetami są i .θ X θ X ~ G m m ( 0,001 , 0,001 ) θ ~ U n i m O r m ( 0 , 1 )
Autor używa niewłaściwego przed ale WinBUGS nie akceptuje niewłaściwych priorytetów.
Przykład
W artykule (strona 226) podano następujące liczby sukcesów zaobserwowanych kozłów wodnych: . Chcę oszacować , liczebność populacji.
Oto jak próbowałem wypracować przykład w WinBUGS ( zaktualizowany po komentarzu @ Stéphane Laurenta):
model {
# Likelihood
for (i in 1:N) {
x[i] ~ dbin(theta, n)
}
# Priors
n ~ dpois(mu)
lambda ~ dgamma(0.001, 0.001)
theta ~ dunif(0, 1)
mu <- lambda/theta
}
# Data
list(x = c(53, 57, 66, 67, 72), N = 5)
# Initial values
list(n = 100, lambda = 100, theta = 0.5)
list(n = 1000, lambda = 1000, theta = 0.8)
list(n = 5000, lambda = 10, theta = 0.2)
Model ma parapecie nie zbiegają się ładnie po 500.000 próbek 20'000 Burn-in próbek. Oto wynik uruchomienia JAGS:
Inference for Bugs model at "jags_model_binomial.txt", fit using jags,
5 chains, each with 5e+05 iterations (first 20000 discarded), n.thin = 5
n.sims = 480000 iterations saved
mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat n.eff
lambda 63.081 5.222 53.135 59.609 62.938 66.385 73.856 1.001 480000
mu 542.917 1040.975 91.322 147.231 231.805 462.539 3484.324 1.018 300
n 542.906 1040.762 95.000 147.000 231.000 462.000 3484.000 1.018 300
theta 0.292 0.185 0.018 0.136 0.272 0.428 0.668 1.018 300
deviance 34.907 1.554 33.633 33.859 34.354 35.376 39.213 1.001 43000
pytania
Najwyraźniej czegoś mi brakuje, ale nie widzę, co dokładnie. Myślę, że moje sformułowanie modelu jest gdzieś błędne. Więc moje pytania to:
- Dlaczego mój model i jego wdrożenie nie działają?
- Jak poprawnie sformułować i wdrożyć model podany przez Raftery (1988)?
Dzięki za pomoc.
źródło
mu=lambda/theta
i wymienićn ~ dpois(lambda)
zn ~ dpois(mu)
Odpowiedzi:
Cóż, skoro masz kod do działania, wygląda na to, że ta odpowiedź jest trochę za późna. Ale już napisałem kod, więc ...
Jeśli chodzi o wartość, jest to ten sam * model, który pasuje(N,θ)
rstan
. Szacuje się to na 11 sekund na moim laptopie konsumenckim, osiągając wyższy efektywny rozmiar próbki dla naszych interesujących parametrów w mniejszej liczbie iteracji.Zauważ, że rzucam
theta
jako 2-simpleks. To tylko dla stabilności numerycznej. Kwota odsetek wynositheta[1]
; oczywiścietheta[2]
jest to zbędna informacja.* Jak widać, tylne podsumowanie jest praktycznie identyczne, a promowanie do rzeczywistej ilości nie wydaje się mieć istotnego wpływu na nasze wnioski.N
Kwantyl 97,5% dla jest o 50% większy dla mojego modelu, ale myślę, że dzieje się tak, ponieważ próbnik Stana jest lepszy w badaniu pełnego zakresu tylnej niż zwykły chód losowy, więc może łatwiej dostać się w reszkę. Jednak mogę się mylić.N
Biorąc wartości wygenerowane ze stanu, używam ich do narysowania późniejszych wartości predykcyjnych ˜ y . Nie powinniśmy się dziwić, że średnia prognoz późniejszych ˜ y jest bardzo zbliżona do średniej z danych przykładowych!N,θ y~ y~
rstan
Poniższy kod może potwierdzić, że nasze wyniki ze stanu mają sens.
rstan
źródło
n
nie można zadeklarować jako liczbę całkowitą i nie znałem obejścia problemu.rstan
wyniki są bardziej poprawne. [1] stats.stackexchange.com/questions/114366/…Oto mój skrypt analizy i wyniki przy użyciu JAGS i R:
Obliczenia trwały około 98 sekund na moim komputerze stacjonarnym.
Wyniki są następujące:
źródło