# Anzahl der Simulationsdurchläufe
<- 10000
R
# Initialisierung eines Vektors, um die
# Werte der Teststatistik aufzunehmen
<- rep(0, R)
teststat
# Parameter
<- 0.05
alpha <- 63
mu0 <- mu0
mu <- 4
sigma <- 100
n
# Schleife
for(r in 1:R){
<- rnorm(n, mean=mu, sd=sigma)
x <- sqrt(n) * (mean(x) - mu0)/sd(x)
teststat[r]
}
15 Bootstrap
Die bisherige Darstellung der Hypothesentests hat möglicherweise den Eindruck erweckt, dass es eine Reihe von vorgefertigten Testverfahren gibt, die man als Nutzer heraussuchen und anwenden kann. Tatsächlich stellt sich in vielen ökonomischen Anwendungen und vor allem in der empirischen akademischen Forschung oft die Frage, wie man eine Hypothese testen kann, für die es noch keinen vorgefertigten Test gibt. Ökonomen sollten also nicht nur vorgefertigte Tests verstehen und anwenden können, sondern auch in der Lage sein, neue Tests zu entwickeln. In diesem Kapitel wird ein sehr flexibler, genereller, computer-gestützter Ansatz für Hypothesentests vorgestellt, bei dem Monte-Carlo-Simulationen eingesetzt werden. Man nennt diese Methode Bootstrap (oder Bootstrapping). Die Bezeichnung bezieht sich auf die englische Redensart “to pull oneself up by one’s bootstraps”. Warum das ein passender Name ist, wird im folgenden deutlich.
15.1 Monte-Carlo-Simulation der kritischen Grenzen
Um die Grundidee des Bootstraps zu beschreiben, betrachten wir ein konkretes Testbeispiel, das in Abschnitt 13.1 bereits ausführlich besprochen wurde, nämlich einen Erwartungswerttest mit den Hypothesen
Die Teststatistik für diesen Test ist
Konkret nehmen wir an, dass
Die Verteilung der Teststatistik unter der Nullhypothese kann man nun durch eine Monte-Carlo-Simulation mit den Methoden aus Abschnitt 8.2 beliebig genau ermitteln. Das prinzipielle Vorgehen ist fast genauso wie in Abschnitt 13.5, in dem wir Fehlerwahrscheinlichkeiten für den Erwartungswerttest durch Monte-Carlo-Simulationen bestimmt haben. Wir interessieren uns jetzt jedoch für das
Der R-Code für die Monte-Carlo-Simulation sieht so aus:
Wenn die Schleife durchgelaufen ist, enthält der Vektor teststat
die R
Realisationen der Teststatistik (unter der Nullhypothese). Die gesuchten Quantile kann man nun mit Hilfe der R-Funktion quantile
bestimmen. Sie sind
quantile(teststat, prob=c(alpha/2, 1-alpha/2))
2.5% 97.5%
-2.022429 1.972603
Diese Werte liegen sehr nah an den theoretisch hergeleiteten kritischen Grenzen, die sich aus der t-Verteilung ergeben, nämlich
qt(c(alpha/2, 1-alpha/2), df=n-1)
[1] -1.984217 1.984217
Wenn die Anzahl der Schleifendurchläufe erhöht wird, dann sind die Unterschiede noch kleiner.
Es ist also mit Hilfe von Monte-Carlo-Simulationen möglich, die kritischen Werte aus der Verteilung der Teststatistik unter der Nullhypothese zu finden, wenn man die wahre Verteilung von
15.2 Stichprobe und Pseudo-Stichproben
Die Idee der Bootstrap-Methode ist eigentlich naheliegend: Da wir die Verteilung von
Ein Problem muss noch gelöst werden: Die Nullhypothese wird von der Original-Stichprobe im allgemeinen nicht exakt erfüllt, möglicherweise wird sie nicht einmal ungefähr erfüllt. Eine Simulation würde also nicht die gesuchte Verteilung der Teststatistik unter Gültigkeit der Nullhypothese ermitteln. Daher muss die Original-Stichprobe so angepasst werden, dass sie die Nullhypothese erfüllt. Dabei soll die Verteilung der Original-Stichprobe möglichst wenig verändert werden.
Für den Erwartungswerttest sieht die Anpassung so aus, dass alle Stichprobenelemente der Original-Stichprobe, also
so dass die Nullhypothese erfüllt ist.
Nun ermittelt man die Verteilung der Teststatistik mit Hilfe von Monte-Carlo-Simulationen unter der Annahme, dass die Populationsverteilung nicht
Für die Bootstrap-Simulation muss in jedem Schleifendurchlauf eine neue Stichprobe aus der Bootstrap-Population sample
ziehen.
<- x - mean(x) + mu0
x0 <- sample(x0, size=n, replace=TRUE) xb
Man zieht also aus der an die Nullhypothese angepassten Original-Stichprobe wie aus einer Urne size=n
). Nach jedem Zug wird das gezogene Element wieder zurückgelegt (replace=TRUE
), so dass ein Stichprobenelement der Original-Stichprobe durchaus mehrfach in der Pseudo-Stichprobe landen kann. Andere Elemente der Original-Stichprobe tauchen dafür in der Pseudo-Stichprobe gar nicht auf.
15.3 Bootstrap-Simulation
Die Verteilung der Teststatistik unter der Nullhypothese wird durch eine Bootstrap-Simulation mit vielen Schleifendurchläufen ermittelt. Die Anzahl der Durchläufe nennen wir
Nachdem die Schleife
Wir gehen nun den R-Code für die Bestimmung der kritischen Grenzen beim Erwartungswerttest Schritt für Schritt durch. Der hypothetische Wert sei x
. Die bootstrap1.csv
unter dem Namen originalstichprobe
abgespeichert.
<- read.csv("../data/bootstrap1.csv")
D head(D)
originalstichprobe
1 58.62
2 56.98
3 63.83
4 61.51
5 64.87
6 54.45
Die Daten werden in dem Vektor x
abgelegt.
<- D$originalstichprobe
x <- length(x) n
Aus der Original-Stichprobe wird die Original-Teststatistik berechnet. Ihr Wert ist
<- 63
mu0 <- sqrt(n)*(mean(x) - mu0)/sd(x)
teststat print(teststat)
[1] -2.348181
Vor dem Start der Bootstrap-Schleife wird ein Vektor teststatb
initialisiert, in dem die Werte der Bootstrap-Teststatistik gesammelt werden können. Die Zahl der Bootstrap-Wiederholungen setzen wir auf
<- 10000
B <- rep(0, B) teststatb
Nun folgt die eigentliche Bootstrap-Routine. In jedem Schleifendurchlauf wird eine Pseudo-Stichprobe xb
aus der an die Nullhypothese angepassten konkreten Stichprobe x0
gezogen. Aus ihr wird die Teststatistik berechnet und als ein Element in dem Vektor teststatb
abgespeichert.
for(b in 1:B){
<- x - mean(x) + mu0
x0 <- sample(x0, size=n, replace=TRUE)
xb <- sqrt(n)*(mean(xb) - mu0)/sd(xb)
teststatb[b]
}
Aus den quantile
-Funktion. Das Signifikanzniveau setzen wir auf 0.05.
<- 0.05
alpha quantile(teststatb,
prob=c(alpha/2, 1-alpha/2))
2.5% 97.5%
-2.104623 1.956397
Der Wert der Original-Teststatistik (-2.35) liegt nicht zwischen diesen beiden Werten, sondern ist kleiner als die untere kritische Grenze. Die Nullhypothese wird also abgelehnt.
Eine Besonderheit des Bootstrapverfahrens besteht darin, dass die kritischen Grenzen speziell für die Original-Stichprobe bestimmt werden. Eine andere Original-Stichprobe führt zu anderen kritischen Grenzen. Natürlich gibt es auch einen Zufallseinfluss durch das Ziehen der Pseudo-Stichproben. Er wird jedoch (aufgrund des Gesetzes der großen Zahl) immer kleiner, wenn man die Zahl der Durchläufe
15.4 Weitere Bootstrap-Beispiele
In den folgenden Abschnitten wird gezeigt, wie einige weitere Hypothesentests mit Hilfe der Bootstrap-Methode durchgeführt werden können. Die einzige Schwierigkeit der Bootstrap-Tests besteht darin, die Original-Stichprobe so anzupassen, dass die Nullhypothese erfüllt wird.
15.4.1 Zwei Erwartungswerte
Wir verwenden die Notation aus Abschnitt 13.4. Es gibt nun zwei Original-Stichproben, nämlich
Es gibt mehrere naheliegende Möglichkeiten: (i) Man verschiebt die
Zuerst werden die Daten der Original-Stichproben eingelesen.
<- read.csv("../data/bootstrap2.csv")
D head(D)
gruppe wert
1 x 63.29
2 x 62.05
3 y 70.37
4 x 57.92
5 x 58.26
6 x 56.54
Nun zieht man aus diesem Dataframe die beiden Stichproben heraus und speichert sie als x
und y
ab.
<- D$wert[D$gruppe == "x"]
x <- D$wert[D$gruppe == "y"] y
Die Stichprobenumfänge sind
<- length(x)
m m
[1] 50
und
<- length(y)
n n
[1] 80
Aus den Original-Stichproben errechnet man den Wert der Teststatistik
<- (mean(x) - mean(y))/sqrt(var(x)/m + var(y)/n)
teststat teststat
[1] -1.461187
Für die Anpassung an die Nullhypothese werden alle Elemente der
<- x - mean(x) + mean(y)
x0 <- y y0
Der weitere R-Code sieht wie folgt aus:
<- 10000
B <- rep(0, B)
teststatb
for(b in 1:B){
<- sample(x0, size=m, replace=TRUE)
xb <- sample(y0, size=n, replace=TRUE)
yb <- (mean(xb) - mean(yb))/sqrt(var(xb)/m + var(yb)/n)
teststatb[b]
}
Die kritischen Grenzen sind auf einem Niveau von 5 Prozent
<- 0.05
alpha quantile(teststatb,
prob=c(alpha/2, 1-alpha/2))
2.5% 97.5%
-2.033125 1.923120
Die Teststatistik (-1.46) liegt zwischen den kritischen Grenzen, daher wird die Nullhypothese nicht abgelehnt. Die Erwartungswerte sind nicht signifikant unterschiedlich.
15.4.2 Varianz-Test
In diesem Abschnitt wird gezeigt, wie man einen Bootstrap-Test in einer neuen Situation anwenden kann. Wir möchten testen, ob die Varianz von
Wie könnte eine Teststatistik aussehen? Die einfachste Lösung ist die Stichprobenvarianz
Der R-Code wird nun vorgestellt. Zuerst wird die Originalstichprobe eingelesen.
<- read.csv("../data/bootstrap3.csv")
D head(D)
originalstichprobe
1 0.4042
2 1.0354
3 0.3947
4 1.5765
5 1.3795
6 0.3585
Die Werte aus dem Dataframe werden in den Vektor x
geschrieben.
<- D$originalstichprobe x
Die Stichprobenlänge ist
<- length(x)
n n
[1] 100
Der Wert der Teststatistik lautet
<- var(x)
teststat teststat
[1] 0.6274378
Es soll getestet werden, ob die Varianz
<- 10000
B <- rep(0, B)
teststatb
for(b in 1:B){
<- x * sqrt(0.5/teststat)
x0 <- sample(x0, size=n, replace=TRUE)
xb <- var(xb)
teststatb[b]
}
Aus den Bootstrap-Teststatistiken werden die kritischen Grenzen bestimmt. Das Signifikanzniveau sei 5 Prozent.
<- 0.05
alpha quantile(teststatb,
prob=c(alpha/2, 1-alpha/2))
2.5% 97.5%
0.2738523 0.7611980
Der Wert der Teststatistik (0.6274) liegt zwischen den kritischen Grenzen. Die Nullhypothese wird also nicht abgelehnt.
15.4.3 Unabhängigkeitstest
Die Bootstrap-Version des Unabhängigkeitstests (Abschnitt 14.2) ist bemerkenswert, weil sie nicht nur bei ausreichend großen Stichproben zuverlässig funktioniert, sondern auch, wenn die Faustregel zur Mindestbesetzung der Tabelleneinträge verletzt ist.
Wir verwenden die Notation aus Abschnitt 14.2. Die Träger von
Wie lässt sich die Verteilung der Teststatistik unter der Nullhypothese, d.h. unter Unabhängigkeit, simulieren? Wenn
Als Beispiel reproduzieren wir das Beispiel aus Abschnitt 14.2: Sind die Tagesrenditen von Siemens und Zalando unabhängig voneinander? Der R-Code sieht so aus:
<- read.csv("../data/tagesrenditen.csv")
renditen <- renditen$siemens
s <- renditen$zalando
z <- length(s)
n
<- c(-Inf, -1.4, -0.6, 0.1, 0.9, 1.6, Inf)
A_breaks <- c(-Inf, -3, -1.5, 0, 1.3, 3.1, Inf) B_breaks
Die Vektoren s
und z
enthalten die Original-Stichproben, also die beobachteten Tagesrenditen von Siemens und Zalando. Die Partitionsgrenzen sind in A_breaks
und B_breaks
gespeichert. Sie bleiben in den Bootstrap-Simulationen unverändert.
In der Bootstrap-Schleife zieht man aus s
und z
unabhängig voneinander die Pseudo-Stichproben sb
und zb
. Aus ihnen bestimmt man den Wert der Teststatistik und speichert ihn ab. Da die Berechnung der Teststatistik recht langsam ist, wird die Bootstrap-Schleife nur 1000 Mal durchlaufen.
<- 1000
B <- rep(0, B)
teststatb
for(b in 1:B){
# Ziehung der Pseudo-Stichproben
<- sample(s, size=n, replace=TRUE)
sb <- sample(z, size=n, replace=TRUE)
zb
# Diskretisierung
<- cut(sb, A_breaks)
sb_d <- cut(zb, B_breaks)
zb_d
# Berechnung der Teststatistik
<- table(sb_d)
Nb_siemens <- table(zb_d)
Nb_zalando <- table(sb_d, zb_d)
Nb <- 0
teststatb[b] for(j in 1:6){
for(k in 1:6){
<- Nb_siemens[j]*Nb_zalando[k]/n
aux <- teststatb[b] + (Nb[j,k]-aux)^2/(aux)
teststatb[b]
}
} }
Gesucht ist nun das obere
<- 0.05
alpha quantile(teststatb, 1-alpha)
95%
36.56356
Der Wert der Teststatistik beträgt 130.9, er liegt also deutlich im kritischen Bereich. Die Nullhypothese wird verworfen. Die Tagesrenditen hängen statistisch signifikant miteinander zusammen.
15.4.4 Anpassungstest
Die Bootstrap-Version des Anpassungstests (Abschnitt 14.3) ist ist auch dann anwendbar, wenn die Faustregel zur Mindestbesetzung der Partitionseinträge nicht erfüllt ist. Wir verwenden die Notation aus Abschnitt 14.3. Die Träger von
wobei sich die theoretischen Wahrscheinlichkeiten
Im Gegensatz zu Pseudo-Stichproben, die aus der Original-Stichprobe gezogen wurden, bietet es sich beim Anpassungstest an, die Pseudo-Stichproben einfach aus der Verteilung zu ziehen, die in der Nullhypothese behauptet wird.
Zur Illustration wiederholen wir ein Beispiel aus Abschnitt 14.3. Die Nullhypothese behauptet, dass die Anzahl der Kunden, die innerhalb einer Viertelstunde bei einer Hotline anrufen, einer Poisson-Verteiung mit dem Parameter
<- 0:9
x_j <- c(6,27,78,81,72,64,33,23,9,7)
n_j
<- length(x_j)
J <- sum(n_j) n
Es gab also 7 Viertelstunden mit jeweils 9 Anrufen, 9 Viertelstunden mit 8 Anrufen etc. Die theoretischen Wahrscheinlichkeiten für die Ausprägungen “0” bis “8” sowie “9 oder mehr” sind
<- c(dpois(0:8, lambda=4),
pi_j 1-ppois(8, lambda=4))
Aus diesen Angaben ergab sich in Abschnitt 14.3 als Original-Teststatistik der Wert 10.3. Wie bestimmt man mit einer Bootstrap-Simulation die kritische Grenze? Der R-Code sieht wie folgt aus:
<- 10000
B <- rep(0, B)
teststatb
for(b in 1:B){
# Ziehung der Pseudo-Stichprobe
<- rpois(n, lambda=4)
xb
# Bestimmung der Häufigkeiten
<- rep(0, J)
Nb for(i in 1:9){
<- sum(xb == (i-1))
Nb[i]
}10] <- sum(xb >= 9)
Nb[
# Berechnung der Teststatistik
<- sum((Nb - n*pi_j)^2/(n*pi_j))
teststatb[b]
}
Zur Erklärung: In der Schleife wird zuerst die Pseudo-Stichprobe aus der hypothetischen Poisson-Verteilung gezogen. Anschließend werden die Häufigkeiten der 10 Ausprägungen “0” bis “8” und “9 oder mehr” ausgezählt und der Wert der Teststatistik berechnet. Aus den
<- 0.05
alpha quantile(teststatb, prob=1-alpha)
95%
17.04704
Da der Wert der Original-Teststatistik mit 10.3 kleiner ist als der kritische Wert, wird die Nullhypothese nicht verworfen. Die Hypothese, dass die Anzahl der Anrufe innerhalb einer Viertelstunde einer Poisson-Verteilung mit dem Parameter
15.5 Fazit
Die Bootstrap-Methode ist ein sehr nützliches Instrument für die statistische Inferenz. Man kann mit Hilfe von Computersimulationen die Verteilung der Teststatistik unter der Nullhypothese simulieren und auf diese Weise den kritischen Bereich bestimmen. Auf welche Weise das genau geschehen soll, hängt zum Teil vom Einzelfall ab. Daher ist es sehr wichtig, dass man ein tiefes Verständnis der statistischen Inferenz entwickelt. Dazu gehört auch, dass man die klassischen Testverfahren gründlich durchdacht hat. Ein reines Anwenden von “Kochrezepten” führt oft in die Irre und bringt auch keinen Spaß.