groesse02

2023

bayes

bayesbox

string

Published

November 8, 2023

Aufgabe

Wir interessieren uns für die typische Körpergröße deutscher Studentis. Hier findet sich dazu ein Datensatz.

Ausgehend von der Annahme, dass sich die Körpergröße normalverteilt (innerhalb eines Geschlechts) suchen wir die Parameter der Normalverteilung, also Mittelwert und Streuung auf Basis unserer Daten. Anders gesagt: Wir suchen die Post-Verteilung für die Körpergröße.

Gehen wir von folgenden Apriori-Wahrscheinlichkeiten für die Parameter der Normalverteilung aus:

Mittelwert: 150cm bis 200 cm, jeder Wert gleich plausibel, alle anderen Werte unmöglich
SD: 1cm bis 20cm, jeder Wert gleich plausibel, alle anderen Werte unmöglich

Ja, das sind ziemlich einfältige Annahmen, aber gut, fangen wir damit an.

Erstellen Sie eine Bayes-Box!

Hinweise:

Untersuchen Sie den angegebenen Parameterbereich in 1cm-Schritten.
Orientieren Sie sich im Übrigen an den allgemeinen Hinweisen des Datenwerks.

Lösung

library(pradadata)  # für den Datensatz `wo_men`
library(prada)  # für bayesbox, alternativ mit `source`
library(tidyverse)
library(ggpubr)

Daten importieren:

data(wo_men)

Alternativ per URL:

wo_men <- read.csv("https://raw.githubusercontent.com/sebastiansauer/pradadata/master/data-raw/wo_men.csv")

Mittelwert in der Stichprobe:

wo_men |> 
  group_by(sex) |> 
  summarise(height_avg = mean(height, na.rm = TRUE),
            height_sd = sd(height, na.rm = TRUE))

sex	height_avg	height_sd
man	183.1111	9.958082
woman	161.3095	42.782380
NA	NaN	NA

Zur Berechnung der Likelihoods betrachten wir nur jeweils ganze Zentimeter der Körpergröße, also z.B. 180cm, 181cm, 182cm. Auf genauere Auflösung (z.B. 180.1cm, 180.01cm, 180.001cm, …) verzichten wir – das wird sonst einfach zu umfangreich. Anders gesagt: Wir diskretisieren wir die stetige Variable height in Stufen von jeweils 1cm, der Einfachheit halber.

Die Wahrscheinlichkeit für das 1cm-Intervall um unserem Stichprobenergebnis herem (182.5cm bis 183.5cm), bei z.B. einem Mittelwert von 180cm und einer SD von 10cm, entspricht dann dieser Differenz:

obere_grenze <- pnorm(q = 183 + 0.5, mean = 180, sd = 10)
untere_grenze <- pnorm(q = 183 - 0.5, mean = 180, sd = 10)

obere_grenze

[1] 0.6368307

untere_grenze

[1] 0.5987063

obere_grenze - untere_grenze

[1] 0.03812433

Visualisieren wir uns kurz dieses Intervall.

library(mosaic)
xpnorm(q = c(182.5, 183.5), mean = 180, sd = 10)

[1] 0.5987063 0.6368307

Als nächstes legen wir die Werte für unsere Bayes-Box fest.

norm_mean <- seq(from = 150, to = 200, by = 1)
norm_sd <- seq(from = 1, to = 20, by = 1)

Jetzt bauen wir unsere Bayes-Box.

Wenn wir die Wahrscheinlichkeiten der Parameter für alle Kombinationen aus 51 Mittelwerten und 20 SD-Werten prüfen wollen, wird die Tabelle ganz schön lang:

anzahl_kombinationen <- length(norm_mean) * length(norm_sd)
anzahl_kombinationen

[1] 1020

Mit expand_grid kann man sich eine Tabelle erstellen lassen, die alle Kombinationen zweier Variablen aufschreibt:

bayes_box <-
  expand_grid(norm_mean, norm_sd)

head(bayes_box)

norm_mean	norm_sd
150	1
150	2
150	3
150	4
150	5
150	6

Das sind unsere Parameterwerte: Jede Kombination eines Mittelwerts und einer Streuung ist eine Hypothese. Insgesamt haben wir also 1020 Parameterwerte.

So, bauen wir die Bayes-Box weiter:

L <- pnorm(183.5, mean = bayes_box$norm_mean, sd = bayes_box$norm_sd)

bayes_box2 <-
  bayes_box |> 
  mutate(hyp = 1:anzahl_kombinationen,
         lik = L,
         post_unstand = hyp * lik,
         post_std = post_unstand / sum(post_unstand))

samples <-
  bayes_box2 |> 
  slice_sample(
    n = 1e4,
    weight_by = post_std,
    replace = TRUE
  )

Und jetzt visualisieren:

samples |> 
  ggplot() +
  aes(x = norm_sd,
      y = norm_mean,
      fill = post_std) +
  geom_tile() +
  scale_fill_viridis_c()

Die Stichproben-Postverteilung erlaubt es auch bequem, die einzelnen Parameter der Post-Verteilung jeweils für sich zu visualisieren:

gghistogram(samples,
            x = "norm_mean",
            bins = 15)

ggdensity(samples,
          x = "norm_mean")

ggdensity(samples,
          x = "norm_sd")

Fazit: Wir haben in dieser Übung die Post-Verteilung des Mittelwerts bestimmt und die Streuung der Körpergrößen.

--- exname: groesse02 expoints: 1 extype: string exsolution: NA categories: - 2023 - bayes - bayesbox - string date: '2023-11-08' slug: groesse02 title: groesse02 ---  ```{r global-knitr-options, include=FALSE, message=FALSE} knitr::opts_chunk$set(fig.pos = 'H', fig.asp = 0.618, fig.width = 4, fig.cap = "", fig.path = "", echo = TRUE, # ECHO IS FALSE!!! message = FALSE, warning = FALSE, fig.show = "hold") ``` # Aufgabe Wir interessieren uns für die typische Körpergröße deutscher Studentis. [Hier](https://raw.githubusercontent.com/sebastiansauer/pradadata/master/data-raw/wo_men.csv) findet sich dazu ein Datensatz. Ausgehend von der Annahme, dass sich die Körpergröße normalverteilt (innerhalb eines Geschlechts) suchen wir die Parameter der Normalverteilung, also Mittelwert und Streuung auf Basis unserer Daten. Anders gesagt: Wir suchen die Post-Verteilung für die Körpergröße. Gehen wir von folgenden Apriori-Wahrscheinlichkeiten für die Parameter der Normalverteilung aus: - Mittelwert: 150cm bis 200 cm, jeder Wert gleich plausibel, alle anderen Werte unmöglich - SD: 1cm bis 20cm, jeder Wert gleich plausibel, alle anderen Werte unmöglich Ja, das sind ziemlich einfältige Annahmen, aber gut, fangen wir damit an. Erstellen Sie eine Bayes-Box! Hinweise: - Untersuchen Sie den angegebenen Parameterbereich in 1cm-Schritten. - Orientieren Sie sich im Übrigen an den [allgemeinen Hinweisen des Datenwerks](https://datenwerk.netlify.app/hinweise). # Lösung ```{r} library(pradadata) # für den Datensatz `wo_men` library(prada) # für bayesbox, alternativ mit `source` library(tidyverse) library(ggpubr) ``` Daten importieren: ```{r} data(wo_men) ``` Alternativ per URL: ```{r} wo_men <- read.csv("https://raw.githubusercontent.com/sebastiansauer/pradadata/master/data-raw/wo_men.csv") ``` Mittelwert in der Stichprobe: ```{r} wo_men |> group_by(sex) |> summarise(height_avg = mean(height, na.rm = TRUE), height_sd = sd(height, na.rm = TRUE)) ``` Zur Berechnung der Likelihoods betrachten wir nur jeweils ganze Zentimeter der Körpergröße, also z.B. 180cm, 181cm, 182cm. Auf genauere Auflösung (z.B. 180.1cm, 180.01cm, 180.001cm, ...) verzichten wir -- das wird sonst einfach zu umfangreich. Anders gesagt: Wir *diskretisieren* wir die stetige Variable `height` in Stufen von jeweils 1cm, der Einfachheit halber. Die Wahrscheinlichkeit für das 1cm-Intervall um unserem Stichprobenergebnis herem (182.5cm bis 183.5cm), bei z.B. einem Mittelwert von 180cm und einer SD von 10cm, entspricht dann dieser Differenz: ```{r} obere_grenze <- pnorm(q = 183 + 0.5, mean = 180, sd = 10) untere_grenze <- pnorm(q = 183 - 0.5, mean = 180, sd = 10) obere_grenze untere_grenze obere_grenze - untere_grenze ``` Visualisieren wir uns kurz dieses Intervall. ```{r} library(mosaic) xpnorm(q = c(182.5, 183.5), mean = 180, sd = 10) ``` Als nächstes legen wir die Werte für unsere Bayes-Box fest. ```{r} norm_mean <- seq(from = 150, to = 200, by = 1) norm_sd <- seq(from = 1, to = 20, by = 1) ``` Jetzt bauen wir unsere Bayes-Box. Wenn wir die Wahrscheinlichkeiten der Parameter für alle Kombinationen aus `r length(norm_mean)` Mittelwerten und `r length(norm_sd)` SD-Werten prüfen wollen, wird die Tabelle ganz schön lang: ```{r} anzahl_kombinationen <- length(norm_mean) * length(norm_sd) anzahl_kombinationen ``` Mit `expand_grid` kann man sich eine Tabelle erstellen lassen, die alle Kombinationen zweier Variablen aufschreibt: ```{r} bayes_box <- expand_grid(norm_mean, norm_sd) head(bayes_box) ``` Das sind unsere Parameterwerte: Jede Kombination eines Mittelwerts und einer Streuung ist eine Hypothese. Insgesamt haben wir also `r nrow(bayes_box)` Parameterwerte. So, bauen wir die Bayes-Box weiter: ```{r} L <- pnorm(183.5, mean = bayes_box$norm_mean, sd = bayes_box$norm_sd) bayes_box2 <- bayes_box |> mutate(hyp = 1:anzahl_kombinationen, lik = L, post_unstand = hyp * lik, post_std = post_unstand / sum(post_unstand)) ``` ```{r} samples <- bayes_box2 |> slice_sample( n = 1e4, weight_by = post_std, replace = TRUE ) ``` Und jetzt visualisieren: ```{r} samples |> ggplot() + aes(x = norm_sd, y = norm_mean, fill = post_std) + geom_tile() + scale_fill_viridis_c() ``` Die Stichproben-Postverteilung erlaubt es auch bequem, die einzelnen Parameter der Post-Verteilung jeweils für sich zu visualisieren: ```{r} gghistogram(samples, x = "norm_mean", bins = 15) ``` ```{r} ggdensity(samples, x = "norm_mean") ``` ```{r} ggdensity(samples, x = "norm_sd") ``` Fazit: Wir haben in dieser Übung die Post-Verteilung des Mittelwerts bestimmt und die Streuung der Körpergrößen.