JAGS: Just Another Gibbs Sampler
Aus WeissWiki
Bayesianische Verfahren gewinnen dank zunehmender Softwareunterstützung an Bedeutung. Wer also nicht in der Lage ist oder keine Lust hat, sich in C oder C++ (oder Fortran oder ...) den eigenen gibbs sampler zu schreiben, der greift häufig auf BUGS beziehungsweise die mit einem GUI versehenen WinBUGS oder OpenBUGS zurück. Vergleichsweise neu ist JAGS (Just Another Gibbs Sampler), das sich eines BUGS Dialektes bedient. JAGS wurde unter anderem mit dem Ziel erstellt, auch unter Linux BUGS-Code verwenden zu können. Darüber hinaus meine ich mich erinnern zu können, dass es schneller als WinBUGS/OpenBUGS ist.
[Bearbeiten] Beispiel: Lineares Modelle mit JAGS und R schätzen
Da R mein Programm der Wahl zur statistischen Datenanalysen ist, liegt es nahe, JAGS aus R heraus nutzen zu wollen. Das heißt insbesondere, (1) Spezifikation des JAGS-Modells, (2) Datenverwaltung und -zuweisung, (3) Steuerung der Stichprobenziehung aus der Posteriorverteilung (erzeugt via MCMC) und (4) Weiterverarbeitung der Ergebnisse, etwa zur Konvergenzdiagnostik (mit Hilfe von coda).
Ich verwende R und JAGS unter MS-Windows.
Noch ein Hinweis: Sämtlicher R Code findet sich auch in folgendem Skript anaJAGSinR.R.
1. Schritt: JAGS herunterladen und installieren
Das Windows binary von JAGS herunterladen und installieren (was dank Installerprogramm kein Problem sein sollte).
2. Schritt: rjags herunterladen und installieren
Die Schnittstelle zwischen JAGS und R bietet das package rjags. Die einfachste Möglichkeit der Installation von rjags stellt die R Funktion install.packages() dar. R starten und in den R Konsole install.packages("rjags") eingeben. Es ist darauf zu achten, dass vorher JAGS installiert wurde (siehe Schritt 1), da rjags nach einer JAGS-Installation sucht.
3. Schritt: JAGS Modell definieren
Nachdem die notwendigen R packages gerladen wurden, kann nun das JAGS Modell definiert werden, dessen Aufbau sich sehr an BUGS anlehnt. Ich verwende hier zur Illustration das 1. Beispiel aus dem JAGS-Manual, das die bayesianische Variante einer linearen Regression ist.
### R packages laden library(rjags) library(BRugs) library(coda)
### Arbeitspfad definieren path <- "d:/tmp/"
### JAGS Modell spezifizieren
modelLM <-function() {
for (i in 1:N) {
Y[i] ~ dnorm(mu[i], tau)
mu[i] <- alpha + beta * (x[i] - x.bar)
}
x.bar <- mean(x)
alpha ~ dnorm(0.0, 0.0001)
beta ~ dnorm(0.0, 0.0001)
sigma <- 1.0/sqrt(tau)
tau ~ dgamma(0.001, 0.001)
}
4. Schritt: Modell speichern, Daten erzeugen und JAGS-Modellobjekt erstellen
JAGS kann das oben definierte Modell nur in Form einer externen Datei einlesen. Daher muss es zunächst in eine Datei geschrieben werden, die anschließend wiederum eingelesen wird. JAGS stellt keine Funktion bereit, um ein in R definiertes JAGS/BUGS-Modell in eine externe Datei zu schreiben. Doch gibt es im R package BRugs die Funktion writeModel, die genau das leistet (eben für BUGS-Modelle). Die Daten für JAGS werden als Liste mydata übergeben. Schließlich wird mit jags.model() das Modell aus der externen Datei wieder eingelesen und es wird der dazugehörige Datensatz definiert.
### JAGS Modell in externe Datei schreiben writeModel(modelLM, con = paste(path,"modelLM.txt",sep = "")) ### Datensatz definieren, 5 Fälle mydata <- list(x = c(1, 2, 3, 4, 5), Y = c(1, 3, 3, 3, 5), N = 5) ### Modell wieder in JAGS einlesen und spezifizieren modelLMc <- jags.model(file = paste(path,"modelLM.txt",sep = ""), data = mydata)
5. Schritt: Stichprobe(n) aus der Posteriorverteilung ziehen
Mit Hilfe der Funktion jags.sample() wird die eigentliche Stichprobenziehung (via MCMC) vorgenommen. Dazu muss der Funktion das JAGS-Modellobjekt übergeben werden, eine Liste von Variablen, für die eintsprechende Stichproben gezogen werden sollen und die Größe der Stichprobe (hier N = 1000).
modelLMs <- jags.samples(model = modelLMc, c("alpha","beta","tau"),n.iter=1000)
Da die Ergebnisse der Simulation anschließend mit coda weiterverarbeitet werden sollen, bietet es sich an, mit Hilfe der Funktion coda.samples() gleich ein coda-konformes Objekt zu erzeugen (also, anstelle von jags.sample() wird coda.sample() genutzt).
modelLMcoda <- coda.samples(model = modelLMc, c("alpha","beta","tau"),n.iter=1000)
6. Schritt: Ergebnisse der Stichprobenziehung mit coda prüfen und Konvergenzdiagnostik
Schließlich lassen sich mit Hilfe der summary-Funktion des coda packages zentrale Lage- und Streuungsstatistiken der Stichprobe ausgeben. Der Aufruf
summary(modelLMcoda)
führt zu folgender Ausgabe in R
> summary(modelLMcoda)
Iterations = 1:1000 Thinning interval = 1 Number of chains = 1 Sample size per chain = 1000
1. Empirical mean and standard deviation for each variable, plus standard error of the mean:
Mean SD Naive SE Time-series SE alpha 3.013 0.5205 0.01646 0.01781 beta 0.782 0.3702 0.01171 0.01126 tau 1.831 1.5254 0.04824 0.06672
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5% alpha 2.052671 2.7451 3.0167 3.2662 4.043 beta 0.009155 0.5963 0.7926 0.9653 1.513 tau 0.153346 0.7364 1.4391 2.4489 6.235
Die Schätzung mit Hilfe des OLS Schätzers lautet wie folgt (Variable x wird zentriert):
> summary(tmp <- lm(Y ~ scale(x, scale = FALSE), data = mydata))
Call: lm(formula = Y ~ scale(x, scale = FALSE), data = mydata)
Residuals: 1 2 3 4 5 -4.00e-01 8.00e-01 4.84e-17 -8.00e-01 4.00e-01
Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 3.0000 0.3266 9.186 0.00273 ** scale(x, scale = FALSE) 0.8000 0.2309 3.464 0.04052 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.7303 on 3 degrees of freedom Multiple R-squared: 0.8, Adjusted R-squared: 0.7333 F-statistic: 12 on 1 and 3 DF, p-value: 0.04052
Noch zeigen sich leichte Unterschiede in den beiden Schätzungen, doch eine Erhöhung der Durchläufe des MCMC Schätzers auf 10.000 Iterationen führt dazu, dass nahezu identische Ergebnisse herauskommen. Abschließend noch eine Grafik zur Konvergenzdiagnostik; sieht nicht berauschend aus. 1000 Iterationen sind also zu wenig.
[Bearbeiten] Tipps & Tricks
[Bearbeiten] Datentransformationen
In Open-/WinBUGS lassen sich Datentransformationen durchführen, etwa die folgende:
for(iin1:N){
z[i] <- sqrt(y[i])
z[i] ~ dnorm(mu,tau)
}
Das ist mit JAGS nicht möglich! sqrt(y[i]) muss also außerhalb von JAGS berechnet werden (siehe auch JAGS-Manual, Kap. 7.0.5).
