Code
library(tidyverse)
library(easystats)
library(ggpubr) # optional6 Die Post befragen
6.1 Lernsteuerung
6.1.1 Position im Modulverlauf
Abbildung 1.1 gibt einen Überblick zum aktuellen Standort im Modulverlauf.
6.1.2 Lernziele
Nach Absolvieren des jeweiligen Kapitels sollen folgende Lernziele erreicht sein.
Sie können …
- die Post-Verteilung anhand einer Stichprobenverteilung auslesen
- Fragen nach Wahrscheinlichkeitsanteilen der Post-Verteilung anhand der Stichprobenverteilung beantworten
- Fragen nach Quantilen anhand der Stichprobenverteilung beantworten
6.1.3 Begleitliteratur
Der Stoff dieses Kapitels orientiert sich an McElreath (2020), Kap. 3.1 und 3.2.
6.1.4 Vorbereitung im Eigenstudium
6.1.5 Benötigte R-Pakete
6.1.6 Begleitvideos
6.2 Zur Erinnerung: Gitterwerte in R berechnen
Berechnen wir mit der Gittermethode (“Bayes-Box”) die Postverteilung für den Globusversuch.
Die Gittermethode ist ein Weg, die Posteriori-Verteilung zu berechnen. Die Posteriori-Verteilung birgt viele nützliche Informationen.
Modell: \(W=6\) Wasser, \(N=9\) Würfen, bei Apriori-Indifferenz gegenüber den Parameterwerten. Und sagen wir \(k=11\) Gitterwerten1, also mit 10 Wasseranteilswerten zwischen 0 und 1.
Übungsaufgabe 6.1 (Welcher Paramterwert hat die höchste Posteriori-Wahrscheinlichkeit?)
n_success <- 6
n_trials <- 9
p_grid <- seq(from = 0, to = 1, by = .1)
L <- dbinom(n_success, size = n_trials, prob = p_grid)
d <-
tibble(p_grid = p_grid,prior = 1) %>%
mutate(likelihood = L) %>%
mutate(unstand_post = (likelihood * prior),
post = unstand_post / sum(unstand_post))- Sequenz von 0 bis 1 mit Schritten der Größe 0.1.
- Likelihood mit 6 Treffern bei 9 Würfen und das Ganze jeweils für alle 11 Parameterwerte
- Dann packen wir alles in eine Tabelle.
Code
library(ggpubr)
ggline(d, x = "p_grid", y = "post") Abb. Abbildung 6.1 zeigt die resultierende Bayes-Box; vor allem ist die Post-Verteilung wichtig.
Voilà, die Post-Verteilung als Tabelle, auch “Bayes-Box” (oder Bayes-Gitter) genannt: s. Tabelle 6.1.
| p_grid | prior | likelihood | unstand_post | post |
|---|---|---|---|---|
| 0.0 | 1 | 0.00 | 0.00 | 0.00 |
| 0.1 | 1 | 0.00 | 0.00 | 0.00 |
| 0.2 | 1 | 0.00 | 0.00 | 0.00 |
| 0.3 | 1 | 0.02 | 0.02 | 0.02 |
| 0.4 | 1 | 0.07 | 0.07 | 0.07 |
| 0.5 | 1 | 0.16 | 0.16 | 0.16 |
| 0.6 | 1 | 0.25 | 0.25 | 0.25 |
| 0.7 | 1 | 0.27 | 0.27 | 0.27 |
| 0.8 | 1 | 0.18 | 0.18 | 0.18 |
| 0.9 | 1 | 0.04 | 0.04 | 0.04 |
| 1.0 | 1 | 0.00 | 0.00 | 0.00 |
6.2.1 Bayes-Box automatisiert
Übrigens kann man die Berechnung der Bayes-Box auch automatisieren, s. Tabelle 6.2 und Listing 6.2, z.B. so:2
Entweder so:
Code
source("https://raw.githubusercontent.com/sebastiansauer/prada/master/R/NAME_bayesbox.R")
bayesbox(hyps = p_grid, priors = 1, liks = L)Mit source importiert man eine R-Skriptdatei. In diesem Fall steht dort der Code für die Funktion bayesbox.
Oder Sie starten das R-Paket, wo die Funktion wohnt:
bayesbox, auch im Paket prada erhältlich
library(prada)
bayesbox(hyps = p_grid, priors = 1, liks = L)bayesbox
- Das Paket
pradasteht nicht im Standard-R-App-Store (“CRAN”), sondern auf Github. Sie könnnen es so installieren:devtools::install_github("sebastiansauer/prada"). - Die Funktion verhält sich wie eine gewöhnliche Bayes-Box: Bei
hypsschreibt man die Hypothesen (bzw. Parameterwerte) auf. Beipriorsdie Priori-Werte und beiliksdie Likelihoods der jeweiligen Hypothesn.
Viele nützliche Fragen (und Antworten) leiten sich ab aus Abb. Abbildung 6.1.
Beispiel 6.1 (Beispiele für Fragen an die Post-Verteilung)
- Mit welcher Wahrscheinlichkeit liegt der Parameter unter einem bestimmten Wert?
- Mit welcher Wahrscheinlichkeit liegt der Parameter zwischen zwei bestimmten Werten?
- Mit 5% Wahrscheinlichkeit liegt der Parameterwert nicht unter welchem Wert?
- Welcher Parameterwert hat die höchste Wahrscheinlichkeit?
- Wie ungewiss ist das Modell über die Parameterwerte?
etc. \(\square\)
Solche Fragen kann man in zwei Gruppen aufteilen:
- Fragen zu Parametern
- Fragen zu Wahrscheinlichkeiten
6.2.2 Bayes-Box für komplexe Modelle
Bisher, für einfache Fragestellungen hat unsere Bayes-Box, das heißt die Gittermethode bestens funktioniert: einfach, robust, formschön3. Allerdings: Funktioniert sie auch bei komplexeren Modellen? Schließlich wollen wir ja auch irgendwann Regressionsmodelle berechnen. Angenommen, wir haben ein Regressionsmodell mit 1 Prädiktor, dann haben wir folgende drei Größen4 zu schätzen: \(\beta_0, \beta_1, \sigma\). Hört sich gar nicht so viel an. Aber Moment, wir müssten dann z.B. die Frage beantworten, wie wahrscheinlich die Daten aposteriori sind, wenn z.B. \(\beta_0 = -3.14\) und \(\beta_1 = 2.71\) und \(\sigma = 0.70\). Demnach müssen wir alle Ausprägungen (“Gitterwerte”) der Variablen multiplizieren. Puh, das wird eine große Zahl. Wenn wir für die drei Größen jeweils 10 Ausprägungen annehmen, was wenig ist, kämen wir \(10\cdot10\cdot10= 1000=10^3\) Kombinationen. Bei 100 Ausprägungen wären es schon \(100^3=10^6\) Kombinationen. Das wäre doch eine recht lange Tabelle.5
Bei einer multiplen Regression mit sagen wir 10 Prädiktoren mit jeweils 100 Ausprägungen rechnet das arme R bis zum jüngsten Tag: \(10^{100}\).
🤖 Bitte tue mir das nicht an!
👨🏫 Schon gut, das können wir R nicht zumuten. Wir brauchen eine andere Lösung!
6.3 Mit Stichproben die Post-Verteilung zusammenfassen
6.3.1 Wir arbeiten jetzt mit Häufigkeit, nicht mit Wahrscheinlichkeit
Kurz gesagt: Komplexere Bayes-Modelle können nicht mehr “einfach mal eben” ausgerechnet werden; die Mathematik wird zu rechenintensiv. Glücklicherweiße gibt es einen Trick, der die Sache nicht nur rechnerisch, sondern auch konzeptionell viel einfacher macht. Dieser Trick lautet: Wir arbeiten nicht mehr mit Wahrscheinlichkeiten, sondern mit Häufigkeiten. Praktischerweise werden wir in Kürze einen R-Golem kennenlernen, der das für uns erledigt. Dieser Golem liefert uns Stichproben aus der Post-Verteilung zurück. Lernen wir jetzt also, wie man mit solchen Stichproben umgeht.
Wichtig
Die Post-Verteilung mit Häufigkeiten (d.h. in Stichprobenform) ist viel einfach zu handhaben als das direkte Arbeiten mit Wahrscheinlichkeiten. Daher sind viele R-Funktionen für Bayes auf Stichproben eingestellt.
Die Bayes-Box-Methode6 ist bei größeren Datensätzen (oder größeren Modellen) zu unpraktisch. In der Praxis werden daher andere, schnellere Verfahren verwendet, sog. Monte-Carlo-Markov-Ketten (MCMC). Wie diese Verfahren funktionieren sind aber nicht mehr Gegenstand dieses Kurses. Wir wenden Sie einfach an, freuen uns und lassen es damit gut sein.7
6.3.2 Häufigkeiten sind einfacher als Wahrscheinlichkeiten
Wie gesagt, typische R-Werkzeuge (“R-Golems”) liefern uns die Post-Verteilung in Stichprobenform zurück. Bevor wir uns aber mit diesen R-Werkzeugen beschäftigen, sollten wir uns vertraut machen mit einer Post-Verteilung in Stichprobenform. Erstellen wir uns also einen Tabelle mit Stichprobendaten aus der Posteriori-Verteilung (Tabelle d), s. Listing 6.3.
samples <-
d %>% # nimmt die Tabelle mit Posteriori-Daten,
slice_sample( # Ziehe daraus eine Stichprobe,
n = 1e4, # mit insgesamt n=10000 Zeilen,
weight_by = post, # Gewichte nach Post-Wskt.,
replace = T) %>% # Ziehe mit Zurücklegen
select(p_grid)Die Wahrscheinlichkeit, einen bestimmten Parameterwert (d.h. aus der Spalte p_grid) aus Tabelle d zu ziehen, ist proportional zur Posteriori-Wahrscheinlichkeit (post) dieses Werts. Ziehen mit Zurücklegen hält die Wahrscheinlichkeiten während des Ziehens konstant. Das Argument weight_by legt die Wahrscheinlichkeit fest, mit der eine Zeile gezogen wird. Wir begnügen uns mit der Spalte mit den Wasseranteilswerten (Parameterwerten), p_grid, die anderen Spalten brauchen wir nicht. Das Ergebnis, Tabelle samples, die aus Stichproben aus der Post-Verteilung besteht, ist (in Auszügen) in Tabelle 6.3 dargestellt.
Wenn Sie jetzt denken: “Warum machen wir das jetzt? Brauchen wir doch gar nicht!” - Dann haben Sie Recht. Künftig werden wir aber, wenn wir mit komplexeren Modellen zu tun haben, nur noch mit Post-Verteilungen auf Stichprobenbasis arbeiten, weil es damit viel einfacher ist.
Tabelle 6.3 zeigt die ersten Zeilen der Stichproben aus der Post-Verteilung.
| p_grid |
|---|
| 0.600 |
| 0.500 |
| 0.700 |
| 0.700 |
| 0.600 |
Hier erstmal die ersten 100 gesampelten Gitterwerte (p_grid):
## [1] 0.6 0.5 0.7 0.7 0.6 0.7 0.8 0.8 0.8 0.8 0.6 0.7 0.7 0.7 0.7 0.5 0.4 0.5
## [19] 0.8 0.7 0.3 0.9 0.7 0.6 0.7 0.6 0.7 0.7 0.4 0.4 0.5 0.8 0.7 0.6 0.7 0.6
## [37] 0.7 0.6 0.8 0.7 0.6 0.6 0.6 0.9 0.9 0.8 0.5 0.4 0.7 0.6 0.8 0.6 0.5 0.7
## [55] 0.8 0.7 0.7 0.6 0.6 0.5 0.7 0.4 0.5 0.9 0.5 0.5 0.9 0.6 0.9 0.6 0.6 0.8
## [73] 0.7 0.5 0.5 0.5 0.5 0.4 0.7 0.7 0.5 0.6 0.8 0.5 0.6 0.7 0.6 0.6 0.7 0.6
## [91] 0.4 0.5 0.6 0.9 0.8 0.5 0.9 0.6 0.8 0.8So sieht unsere “Stichproben-Bayesbox” als Balkendiagramm aus, s. Abbildung 6.2.
Code
samples |>
count(p_grid) |>
ggbarplot(x = "p_grid", y = "n") 
Aus Abbildung 6.2 können wir einfach auslesen, wie wahrscheinlich gewisse Parameterwerte sind. So sehen wir, dass das Modell Parameterwerte (Wasseranteil, \(\pi\)) zwischen ca. 50% und 70% für am wahrscheinlichsten hält. Aber auch kleine Anteile wie 25% sind nicht auszuschließen (auf Basis der Daten und der Modellannahmen).
Vergleichen Sie Abbildung 6.2 mit Abbildung 5.12: beide sind sehr ähnlich! Das Stichprobenziehen (Abbildung 6.2) nähert sich recht gut an die exakte Berechnung an (Abbildung 5.12).
Es ist hilfreich, sich die Stichproben-Posterior-Verteilung einmal anzuschauen, um ein Gefühl dafür zu bekommen, wie die Post-Verteilung aussieht.
6.3.3 Visualisierung der Stichprobendaten mit \(k=100\) Gitterwerten
\(k=10\) Gitterwerte ist ein grobes Raster. Drehen wir mal die Auflösung auf \(k=100\) Gitterwerte (Ausprägungen) nach oben.
Code
k <- 100
n_success <- 6
n_trials <- 9
d_k100 <-
tibble(p_grid = seq(from = 0,
to = 1,
length.out = k),
prior = 1) %>%
mutate(likelihood = dbinom(n_success,
size = n_trials,
prob = p_grid)) %>%
mutate(unstand_post = (likelihood * prior),
post = unstand_post / sum(unstand_post))- 1
- \(k=100\) Gitterwerte
- 2
- 6 Treffer (Wasser)
- 3
- 9 Versuche
- 4
- Bayes-Box anlegen mit 100 Zeilen, d.h. 100 Parameterwerten
- 5
- Apriori indifferent: Alle Hypothesen haben die gleiche Apriori-Plausibilität
- 6
- Die Likelihood ist binomialverteilt.
- 7
- Post-Verteilung berechnen wie gewohnt.
\(d_k100\) ist eine Bayes-Box mit \(W=6, N=9, k=100\).
Und daraus ziehen wir uns \(n=1000\) Stichproben:
Code
samples_k100 <-
d_k100 %>% # nimmt die Tabelle mit Posteriori-Daten,
slice_sample( # Ziehe daraus eine Stichprobe,
n = 1000, # mit insgesamt n=1000 Elementen,
weight_by = post, # Gewichte nach Spalte mit Post-Wskt.,
replace = T) # Ziehe mit ZurücklegenAbbildung 6.3 zeigt sowohl die exakte Post-Verteilung als auch die Post-Verteilung auf Basis von Stichproben. Im mittleren Teildiagramm sind die Stichproben einzeln als Kreis dargestellt. Im rechten Teildiagramm sind die gleichen Daten als Dichtediagramm dargestellt. In allen Fällen erkennt man gut die zentrale Tendenz: ein Wasseranteil von 70% scheint der “typische” Wert des Modells zu sein. Außerdem erkennt man, dass das Modell durchaus einige Streuung in der Schätzung des Wasseranteils bereithält. Das Modell ist sich nicht sehr sicher, könnte man sagen.
Die Stichprobendaten nähern sich der “echten” Posteriori-Verteilung an: Die Stichproben-Post-Verteilung hat jetzt “glattere” Ränder.
Hinweis
Mehr Stichproben und mehr Gitterwerte glätten die Verteilung.
Jetzt die Post-Verteilung noch mal mit mehr Stichproben: \(n=10^6\) Stichproben bei \(k=100\) Gitterwerten aus der Posteriori-Verteilung, s. Abbildung 6.4.
6.4 Die Post-Verteilung befragen
So, jetzt befragen wir die Post-Verteilung.
📺 Die Post-Verteilung auslesen
Wichtig
Die Post-Verteilung ist das zentrale Ergebnis einer Bayes-Analyse. Wir können viele nützliche Fragen an sie stellen.
Es gibt zwei Arten von Fragen:
- nach Wahrscheinlichkeiten (p)
- nach Parameterwerten (Quantilen, q)
Der Unterschied zwischen beiden Arten von Fragen ist in Abbildung 6.5 schematisch illustriert.
Im linken Teildiagramm von Abbildung 6.5 fragen wir: “Wie wahrscheinlich ist ein Wasseranteil von höchstens 80%?”. Im rechten Teildiagramm fragen wir: “Welcher Wasseranteil wird mit einer Wahrscheinlichkeit von 78% nicht überschritten?”.
6.4.1 Fragen nach Wahrscheinlichkeiten
Sagen wir, dass sei unsere Forschungsfrage: Wie groß ist die Wahrscheinlichkeit, dass der Wasseranteil unter 50% liegt? Um diese Frage zu beantworten, zählen wir einfach, wie viele Stichproben die Bedingung erfüllen, und summieren die Wahrscheinlichkeiten dieser Stichproben. Wir zählen (count) also die Stichproben, die sich für einen Wasseranteil (p_grid) von weniger als 50% aussprechen.
Code
samples %>%
count(p_grid < .5) Da wir insgesamt 10000 (1e4) Stichproben gezogen haben, können wir noch durch diese Zahl teilen, um einen Anteil zu bekommen. Dieser Anteil ist die Antwort auf die Forschungsfrage: Wie Wahrscheinlichkeit (laut Modell) für einen Wasseranteil kleiner als 50%.
Beispiel 6.2 (Was macht die Funktion count?) Der Befehl count macht Folgendes: Er gruppiert die Stichprobe nach dem Prüfkriterium, Wasseranteil höchstens 50%. Dann zählt er in jeder der beiden Teiltabelle die Zeilen und liefert diese zwei Zahlen dann zurück. \(\square\)
Wir zählen wie oft der Wasseranteil weniger als 50% beträgt.
Code
samples %>%
count(p_grid < .5) Natürlich gibt es verschiedene Wege, die gleiche Frage zu beantworten.
Code
d %>%
filter(p_grid < .5) %>%
summarise(sum = sum(post))Beispiel 6.3 (Wasseranteil zwischen 50 und 75%?) Noch eine Forschungsfrage: Mit welcher Wahrscheinlichkeit liegt der Parameter (Wasseranteil) zwischen 0.5 und 0.75?
Wir zählen die Stichproben, die diesen Kriterien entsprechen.
Code
samples %>%
count(p_grid > .5 & p_grid < .75)🤖 Ich würde empfehlen, die Anzahl noch in Anteile umzurechnen. Die kann man dann als Wahrscheinlichkeiten auffassen.
👨🏫 Das wollte ich auch gerade sagen…
Code
samples %>%
count(p_grid > .5 & p_grid < .75) %>%
summarise(Anteil = n / 1e4,
Prozent = 100 * n / 1e4) # In ProzentAnteile von count() könnte man, wenn man möchte, auch filter() verwenden.
Code
samples %>%
filter(p_grid > .5 & p_grid < .75) %>%
summarise(sum = n() / 1e4,
anteil = 100 * n() / 1e4) # In ProzentFertig 😄 \(\square\)
Beispiel 6.4 (Wasseranteil zwischen 90& und 100%?) Noch ein Beispiel für eine Forschungsfrage: Mit welcher Wahrscheinlichkeit liegt der Parameter zwischen 0.9 und 1?
Code
samples %>%
count(p_grid >= .9 & p_grid <= 1) %>%
summarise(prop = 100 * n() / 1e4) # prop wie "proportion", AnteilLaut unserem Modell ist es also sehr unwahrscheinlich, dass der Wasseranteil der Erde mind. 90% beträgt. \(\square\)
Übungsaufgabe 6.2 (Wasseranteil höchstens 50%?)
👩🔬 Mit welcher Wahrscheinlichkeit ist der Planet höchstens zur Hälfte mit Wasser bedeckt?
Wir können auch fragen, welcher Parameterwert am wahrscheinlichsten ist; dieser Wert entspricht dem “Gipfel” des Berges, s. Abbildung 6.4.
Für unsere Stichproben-Postverteilung, samples, s. Abbildung 6.2, lässt sich der Modus so berechnen:
Code
map_estimate(samples$p_grid) Dabei steht map für Maximum Aposteriori, also das Maximum der Post-Verteilung.
Übungsaufgabe 6.3 Bei der Gelegenheit könnten wir folgende, ähnliche Fragen stellen:
- Was ist der mittlere Schätzwert (Mittelwert) zum Wasseranteil laut Post-Verteilung?
- Was ist der mediane Schätzwert (Median)?
6.4.2 Fragen nach Parameterwerten
Wichtig
Schätzbereiche von Parameterwerten nennt man auch Konfidenz- oder Vertrauensintervall8.
Welcher Parameterwert wird mit 90% Wahrscheinlichkeit nicht überschritten, laut unserem Modell? (Gesucht sind also die unteren 90% der Posteriori-Wahrscheinlichkeit) Wir möchten also ziemlich sicher, was die Obergrenze an Wasser auf diesem Planeten ist9.
Code
samples %>%
summarise(quantil90 = quantile(p_grid, p = .9))Laut unserem Modell können wir zu 90% sicher sein, dass der Wasseranteil kleiner ist als ca. 78%.
Es hilft vielleicht, sich die Post-Verteilung noch einmal vor Augen zu führen, s. Abbildung 6.4.
Was ist das mittlere Intervall, das mit 90% Wahrscheinlichkeit den Parameterwert enthält, laut dem Modell?
Dafür “schneiden” wir links und rechts die 5% der Stichproben mit den extremsten Werten ab und schauen, bei welchen Parameterwerten wir als Grenzwerte landen:
Code
samples %>%
summarise(
quant_05 = quantile(p_grid, 0.05),
quant_95 = quantile(p_grid, 0.95))Solche Fragen lassen sich also mit Hilfe von Quantilen beantworten.
6.4.3 Zur Erinnerung: Quantile
Beispiel: Wie groß sind die Studentis (Quelle des Datensatzes)?
Das Quantil von z.B. 25% zeigt die Körpergröße der 25% kleinsten Studentis an, analog für 50%, 75%, in Inches bzw. in Zentimetern10.
Code
# speed_gender_height <- read.csv("https://raw.githubusercontent.com/rpruim/OpenIntro/master/data/speed_gender_height.csv")
data("speed_gender_height", package = "openintro")
height_summary <-
speed_gender_height %>%
mutate(height_cm = height*2.54) %>%
select(height_inch = height, height_cm) %>%
drop_na() %>%
pivot_longer(everything(), names_to = "Einheit", values_to = "Messwert") %>%
group_by(Einheit) %>%
summarise(q25 = quantile(Messwert, prob = .25),
q50 = quantile(Messwert, prob = .5),
q75 = quantile(Messwert, prob = .75))
height_summary- 1
- Daten importieren
- 2
- Inch in Zentimeter umrechnen
- 3
- Zeilen mit fehlenden Werten löschen
- 4
- In die Langform pivotieren
- 5
- Gruppieren nach Einheit (Inch, Zentimeter)
- 6
- Quantile berechnen (Q1, Q2, Q3)
Das 25%-Quantil nennt man auch 1. Quartil; das 50%-Quantil (Median) auch 2. Quartil und das 75%-Quantil auch 3. Quartil.
Abbildung 6.6 visualisiert die Quantile und die Häufigkeitsverteilung.
Übungsaufgabe 6.4 (Welcher Parameterwert ist der wahrscheinlich größte?) Übersetzen wir “wahrscheinlich” größte in “mit einer Wahrscheinlichkeit von 99% gibt es keinen größeren”.
Übungsaufgabe 6.5 (Welcher Parameterwert ist der wahrscheinlich kleinste?) Übersetzen wir “wahrscheinlich” kleinste in “mit einer Wahrscheinlichkeit von 99% gibt es keinen kleineren”.
Übungsaufgabe 6.6 (Welcher Parameterwert ist der “vermutlich” kleinste?) In der “wirklichen” Welt sind Aussagen nicht immer präzise. Sagen wir, die Chefin der Weltraumbehörde hat in einem Presse-Statement von der “vermutlichen Untergrenze” hinsichtlich des Wasseranteils gesprochen.
Übersetzen wir “vermutlich” kleinste in “mit einer Wahrscheinlichkeit von 90% gibt es keinen kleineren”.
6.4.4 Den Quantilen unter die Motorhaube geschaut
Den R-Befehl quantile() kann man sich, wenn man will, einfach nachbauen und entmystifizieren.
Angenommen, wir wollen wissen, welcher Wasseranteil mit 90% Wahrscheinlichkeit nicht überschritten wird. Das können wir mit im Datensatz samples so erreichen.
- Sortiere die Stichproben aufsteigend.
- Schneide die oberen 10% (von 10000) ab (entferne sie).
- Schaue, was der größte verbleibende Wert ist.
Code
samples %>%
arrange(p_grid) %>% # sortiere
slice_head(n = 9000) %>% # nur die ersten 90%
summarise(p90 = max(p_grid))Das (annähernd) gleiche Ergebnis liefert quantile():
Code
samples %>%
summarise(q90 = quantile(p_grid, .9))6.4.5 Visualisierung der Intervalle
Definition 6.1 (Perzentilintervall (PI)) Intervalle (Bereiche), die die “abzuschneidende” Wahrscheinlichkeitsmasse hälftig auf die beiden Ränder aufteilen, nennen wir Perzentilintervalle oder (synonym) Equal-Tails-Intervalle (ETI), s. Abb. Abbildung 6.7, rechtes Teildiagramm.11 \(\square\)
Das 10%, 20%, … 100%-Quantil12 (auf Basis von samples_k100) sind in Abbildung 6.8 illustriert.
samples_k100
6.5 Schiefe Posteriori-Verteilungen sind möglich
Noch einmal zum Globusversuch: Gehen wir von 3 Würfen mit 3 Mal Wasser (Treffer) aus; auf welche Wasseranteile (Parameterwerte) werden wir jetzt schließen?
Vermutlich ziemlich hohe.
Erstellen wir uns dazu mal eine Post-Verteilung (3 Treffer, 3 Würfe), s. Listing 6.4.
d_33 <-
tibble(p_grid = seq(0,1, by =.01),
prior = 1) %>%
mutate(likelihood = dbinom(3, size = 3, prob = p_grid)) %>%
mutate(unstand_post = likelihood * prior) %>%
mutate(post_33 = unstand_post / sum(unstand_post))
samples_33 <-
d_33 %>%
slice_sample(n = 1e6,
weight_by = post_33,
replace = T)So sehen die ersten paar Zeilen der Post-Verteilung, samples_33, aus.
| p_grid | prior | likelihood | unstand_post |
|---|---|---|---|
| 0.79 | 1 | 0.49 | 0.49 |
| 0.77 | 1 | 0.46 | 0.46 |
| 1.00 | 1 | 1.00 | 1.00 |
| 0.92 | 1 | 0.78 | 0.78 |
| 0.63 | 1 | 0.25 | 0.25 |
| 0.62 | 1 | 0.24 | 0.24 |
Mit dieser “schiefen” Post-Verteilung können wir gut die Auswirkungen auf das Perzentil- und das Höchste-Dichte-Intervall anschauen.
6.5.1 Perzentilintervall
Hier z.B. ein 50%-Perzentilintervall, s. Abb. Abbildung 6.9.
Ein Perzentilintervall (ETI, PI) kann, wenn es dumm läuft, den wahrscheinlichsten Parameterwert nicht enthalten, diesen Wert also plausiblen Wert also zurückweisen. Das ist nicht so toll.
Ein Highest-Density-Intervall (HDI13) ist schmäler als der Perzintilintervall und enthält immer den wahrscheinlichsten Parameterwert.
Die Grenzwerte dieses ETI (oder jedes beliebig breiten) kann man sich z.B. mit dem Befehl eti ausgeben lassen.
Code
samples_33 %>%
select(p_grid) %>%
eti(ci = .5) # Paket `easystats`Der wahrscheinlichste Parameterwert (1) ist nicht im Intervall enthalten. Das ist ein Nachteil der ETI.
6.5.2 Intervalle höchster Dichte
Definition 6.2 Intervalle höchster Dichte (Highest density Intervals, HDI oder HDPI) sind definiert als das schmälste Intervall, das den gesuchten Parameter enthält (in Bezug auf ein gegebenes Modell).
Der wahrscheinlichste Parameterwert (\(1\)) ist im Intervall enthalten, was Sinn macht. Bei einem HDI sind die abgeschnitten Ränder nicht mehr gleich groß, im Sinne von enthalten nicht (zwangsläufig) die gleiche Wahrscheinlichkeitsmasse. Bei PI ist die Wahrscheinlichkeitsmasse in diesen Rändern hingegen gleich groß.
Je symmetrischer die Verteilung, desto näher liegen die Punktschätzer aneinander (und umgekehrt), s. Abb. Abbildung 6.10.
Mit dem Befehl hdi kann man sich die Grenzwerte eines HDI, z.B. eines 50%-HDI, ausgeben lassen, s. Tabelle 6.4.
Code
samples %>%
select(p_grid) %>%
hdi(ci = .5) # aus dem Paket `{easystats}`Das Modell ist sich also zu 50% sicher, dass der gesuchte Parameter (der Wasseranteil der Erdoberfläche) sich im von ca. .67 bis .78 befindet (auf Basis eines HDI).
6.6 Fazit
6.6.1 Intervalle höchster Dichte vs. Perzentilintervalle
- Bei symmetrischer Posteriori-Verteilung sind beide Intervalle ähnlich
- Perzentilintervalle sind verbreiteter
- Intervalle höchster Dichte (Highest Density Interval, HDI) sind bei schiefen Post-Verteilungen zu bevorzugen
- Intervalle höchster Dichte sind die schmalsten Intervalle für eine gegebene Wahrscheinlichkeitsmasse
6.6.2 Zusammenfassung
Fassen wir zentrale Punkte an einem Beispiel zusammen.
Im Globusversuch, Datensatz samples, s. Listing 6.3. Sagen wir, wir haben 6 Treffer bei 9 Würfen erzielt.
Lageparameter: Welchen mittleren Wasseranteil kann man erwarten?
Code
samples %>%
summarise(
mean = mean(p_grid),
median = median(p_grid)) Streuungsparameter: Wie unsicher sind wir in der Schätzung des Wasseranteils?
Code
samples %>%
summarise(
p_sd = sd(p_grid),
p_iqr = IQR(p_grid),
p_mad = mad(p_grid)) # Mean Absolute Deviation, Mittlerer AbsolutfehlerAnstelle der Streuungsparameter ist es aber üblicher, ein HDI oder PI anzugeben.
Wichtig
Alles Wasser oder was? Im Beispiel dieses Kapitels haben wir unser gefragt, was wohl der Wasseranteil auf dem Planeten Erde ist. Halten Sie sich klar vor Augen: Der Wasseranteil ist ein Beispiel für einen Parameter, einer unbekannten Größes eines Modells.
6.7 Aufgaben
6.7.1 Papier-und-Bleistift-Aufgaben
6.7.2 Aufgaben, bei denen man einen Computer benötigt
6.8 —
Die Anzahl der Gitterwerte ist nicht Teil des Modells, streng genommen; die Anzahl der Gitterwerte entscheiden nur über die Genauigkeit der Post-Verteilung.↩︎
Alternativ kann man auf die Funktion über das R-Paket {prada} zugreifen.↩︎
naja, nicht unbedingt formschön, aber mir fiel kein dritter Vorzug ein.↩︎
Modellparameter genannt↩︎
Vorsicht beim Ausdrucken.↩︎
auch Grid-Methode genannt↩︎
Eine gute Einführung in die Hintergründe findet sich bei McElreath (2020).↩︎
Tatsächlich gibt es eine Vielzahl an Begriffen, die in der Literatur nicht immer konsistent verwendet werden, etwa Kompatibilitätsintervall, Ungewissheitsintervall, Passungsbereich.↩︎
Vielleicht damit es genug Berge zum Schifahren gibt.↩︎
1 Inch entspricht 2.54cm↩︎
Hier auf Basis der Post-Verteilung
samples.↩︎d.h. die Dezile↩︎
Auch als Highest Hensity Posterior Interval (HDPI) bezeichnet.↩︎