6 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

library(tidyverse)
library(easystats)
library(ggpubr)  # optional

6.1.6 Begleitvideos

Playlist QM2

6.2 Mit Stichproben die Post-Verteilung zusammenfassen

6.2.1 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\) Gitterwerten¹, also mit 10 Wasseranteilswerten zwischen 0 und 1.

Übungsaufgabe 6.1 (Welcher Paramterwert hat die höchste Posteriori-Wahrscheinlichkeit?)

Lösung

Aufgrund der Apriori-Indifferenz entsprechen die Posteriori-Wahrscheinlichkeiten den Likelihoods. Die höchste Wahrscheinlichkeit (d. h. Likelihood) hat derjenige Parameterwert, zu dem die Daten am besten passen, und das ist 6/9 = 2/3, d. Listing 6.1. \(\square\)

Listing 6.1: Bayes-Box mit 6 Wasser bei 9 Versuchen

n <- 11
n_success <- 6          
n_trials  <- 9
p_grid <- seq(from = 0, to = 1, length.out = n)
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 der Länge n=11 Schritte, das gibt uns die Parameterwerte
Likelihood mit 6 Treffern bei 9 Würfen und das Ganze jeweils für alle 11 Parameterwerte
Dann packen wir alles in eine Tabelle.

Abb. Abbildung 6.1 zeigt die resultierende Bayes-Box; vor allem ist die Post-Verteilung wichtig.

library(ggpubr)
ggline(d, x = "p_grid", y = "post")

Abbildung 6.1: Die Postverteilung für W=6, N=9, k=10

Voilà, die Post-Verteilung als Tabelle, auch “Bayes-Box” (oder Bayes-Gitter) genannt: s. Tabelle 6.1.

Tabelle 6.1: Postverteilung mit der Gittermethode berechnet

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.2 Bayes-Box automatisiert

Übrigens kann man die Berechnung der Bayes-Box auch automatisieren, s. Tabelle 6.2 und Listing 6.2, z. B. so:²

Entweder so:

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:

Listing 6.2: Funktion bayesbox, auch im Paket prada erhältlich

library(prada)
bayesbox(hyps = p_grid, priors = 1, liks = L)

Tabelle 6.2: Eine Bayes-Box ‘automatisiert’ erstellt, mit Hilfe der Funktion bayesbox

Das Paket prada steht 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 hyps schreibt man die Hypothesen (bzw. Parameterwerte) auf. Bei priors die Priori-Werte und bei liks die 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.3 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ön³. 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ößen⁴ 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.⁵

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.2.4 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 umfangreich bzw. zu komplex. 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-Methode⁶ ist bei größeren Datensätzen (oder größeren Modellen) zu unpraktisch. In der Praxis werden d. h.r 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⁷

6.2.5 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.

Listing 6.3: Wir stellen eine Tabelle mit Stichproben aus der Post-Verteilung

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.

Tabelle 6.3: Stichproben-Post-Verteilung

p_grid
0.500
0.900
0.600
0.700
0.500

Hier erstmal die ersten 100 gesampelten Gitterwerte (p_grid):

##   [1] 0.5 0.9 0.6 0.7 0.5 0.7 0.9 0.6 0.5 0.8 0.4 0.6 0.7 0.4 0.8 0.6 0.6 0.5
##  [19] 0.6 0.6 0.8 0.7 0.5 0.6 0.5 0.9 0.4 0.7 0.8 0.7 0.5 0.7 0.8 0.6 0.6 0.7
##  [37] 0.7 0.7 0.7 0.8 0.6 0.6 0.7 0.7 0.8 0.5 0.6 0.4 0.6 0.8 0.7 0.7 0.8 0.7
##  [55] 0.6 0.7 0.4 0.7 0.6 0.7 0.8 0.8 0.7 0.8 0.6 0.5 0.6 0.7 0.6 0.6 0.5 0.7
##  [73] 0.4 0.8 0.6 0.6 0.6 0.7 0.7 0.8 0.7 0.4 0.6 0.6 0.7 0.4 0.5 0.7 0.8 0.5
##  [91] 0.6 0.5 0.6 0.5 0.7 0.8 0.7 0.9 0.7 0.8

So sieht unsere “Stichproben-Bayesbox” als Balkendiagramm aus, s. Abbildung 6.2.

Syntax
Output

samples |> 
  count(p_grid) |> 
  ggbarplot(x = "p_grid", y = "n")

Abbildung 6.2: Stichprobenverteilung auf Basis von Stichproben

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.10: beide sind sehr ähnlich! Das Stichprobenziehen (Abbildung 6.2) nähert sich recht gut an die exakte Berechnung an (Abbildung 5.10).

6.2.6 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.

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:

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ücklegen

Abbildung 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.

Abbildung 6.3: Post-Verteilung mit 100 Gitterwerten, exakt vs. auf Basis von Stichproben

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.

Abbildung 6.4: Post-Verteilung mit vielen Stichproben und vielen Parameterwerten (Gitterwerten): schön ‘glatt’. Mittelwert (MW), Modus und Median (Md) liegen eng nebeneinander, da die Verteilung recht symmetrisch ist.

6.3 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.

Abbildung 6.5: Fragen nach p vs. Fragen nach q

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.3.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:

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:

samples %>%
  count(p_grid < .5)

Man könnte also alternativ auch schreiben:

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:

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…

samples %>% 
  count(p_grid > .5 & p_grid < .75) %>% 
  summarise(Anteil = n / 1e4,
            Prozent = 100 * n / 1e4)  # In Prozent

Anteile von count() könnte man, wenn man möchte, auch filter() verwenden:

samples %>% 
  filter(p_grid > .5 & p_grid < .75) %>% 
  summarise(sum     =       n() / 1e4,
            anteil = 100 * n() / 1e4)  # In Prozent

Fertig 😄 \(\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?

Syntax
Output

samples %>% 
  count(p_grid >= .9 & p_grid <= 1) %>% 
  summarise(prop = 100 * n() / 1e4)  # prop wie "proportion", Anteil

Laut 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?

Lösung

samples %>% count(p_grid <= .5)

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:

map_estimate(samples$p_grid)  
## MAP Estimate: 0.70

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)?

Lösung

samples %>% 
  summarise(mean(p_grid),
            median(p_grid))

6.3.2 Fragen nach Parameterwerten

Wichtig

Schätzbereiche von Parameterwerten nennt man auch Konfidenz- oder Vertrauensintervall⁸.

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 ist⁹.

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:

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.3.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¹⁰:

# 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.

Abbildung 6.6: Größenverteilung von 1325 amerikanischen Studentis

Ü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”.

Lösung

samples |> 
  summarise(quant99 = quantile(p_grid, p = .99))

Mit einer Wahrscheinlichkeit von 99% ist der höchste zu erwartende Wasseranteil 0.9.

Übungsaufgabe 6.5 (Welcher Parameterwert ist der wahrscheinlich kleinste?) Übersetzen wir “wahrscheinlich” kleinste in “mit einer Wahrscheinlichkeit von 99% gibt es keinen kleineren”.

Lösung

samples |> 
  summarise(wahrscheinlich_kleinste = 
              quantile(p_grid, p = .01))

Mit einer Wahrscheinlichkeit von 99% ist der kleinste zu erwartende Wasseranteil 0.3 – immer auf Basis unserer beobachteten Daten und der Vorannahmen.

Ü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”.

Lösung

samples |> 
  summarise(wahrscheinlich_kleinste = 
              quantile(p_grid, p = .1))

Mit einer Wahrscheinlichkeit von 90% ist der kleinste zu erwartende Wasseranteil 0.5 – immer auf Basis unserer beobachteten Daten und der Vorannahmen.

6.3.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.

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():

samples %>% 
  summarise(q90 = quantile(p_grid, .9))

6.3.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.¹¹ \(\square\)

(a) Intervall der Post-Verteilung mit den unteren 80% der Wahrscheinlichkeit

Das 10%, 20%, … 100%-Quantil¹² (auf Basis von samples_k100) sind in Abbildung 6.8 illustriert.

Abbildung 6.8: Quantile in 10%-Schritenn durch die Verteilung von `samples_k100`

6.4 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:

Listing 6.4: Schiefe Post-Verteilung in einer Bayes-Box

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 = 1e4, 
                 weight_by = post_33, 
                 replace = T)

So sehen die ersten paar Zeilen der Post-Verteilung, samples_33, aus.

p_grid	prior	likelihood	unstand_post
1.00	1	1.00	1.00
0.87	1	0.66	0.66
0.61	1	0.23	0.23
0.66	1	0.29	0.29
0.87	1	0.66	0.66
0.76	1	0.44	0.44

Mit dieser “schiefen” Post-Verteilung können wir gut die Auswirkungen auf das Perzentil- und das Höchste-Dichte-Intervall anschauen.

6.4.1 Perzentil-Intervall

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 (HDI¹³) 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. so ausgeben lassen:

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.4.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.

Abbildung 6.10: Visualisierung der Punktschätzer bei einer schiefen Post-Verteilung

So kann man sich die Grenzwerte eines 50%-HDI ausgeben lassen, s. Tabelle 6.4.

samples %>% 
  select(p_grid) %>% 
  hdi(ci = .5)  # aus dem Paket `{easystats}`

Tabelle 6.4: 50%-HDI für unser Globusmodell

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.5 Fazit

6.5.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.5.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?

samples %>% 
  summarise(
    mean   = mean(p_grid),
    median = median(p_grid))

Streuungsparameter: Wie unsicher sind wir in der Schätzung des Wasseranteils?

samples %>% 
  summarise(
    p_sd   = sd(p_grid),
    p_iqr = IQR(p_grid),
    p_mad = mad(p_grid))  # Mean Absolute Deviation, Mittlerer Absolutfehler

Anstelle 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.6 Aufgaben

6.7 —

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.↩︎