###
# code 9
# 240604
###

if (!require("pacman")) install.packages("pacman"); library("pacman")
p_load(data.table)
p_load(ggplot2)
p_load(infer)

set.seed(240604)

# Initialisiere data.table
wähler <- data.table()

# Generiere Personen
wähler[, person_id := 1:100000]

# Generiere Alter mit einer realistischeren Verteilung
wähler[, alter := sample(18:80, size = 100000, replace = T)]

# Generiere Einkommen mit einer lognormalen Verteilung
wähler[, einkommen := round(rlnorm(.N, meanlog = 10, sdlog = 0.5))]

# Generiere Berufstypen
wähler[, beruf := sample(c('Studierende/r', 'Arbeiter/in', 'Angestellte/r', 'Selbständig'), size = .N, replace = TRUE, prob = c(0.1, 0.3, 0.3, 0.3))]

# Generiere Familienstand
wähler[, verheiratet := "Nein"]
wähler[alter > 25, verheiratet := sample(c('Ja', 'Nein'), size = .N, replace = TRUE)]

# Generiere Anzahl der Kinder
wähler[, kinder := 0]
wähler[alter > 30, kinder := sample(0:3, size = .N, replace = TRUE)]

# Generiere Verkehrsmittel
wähler[, verkehrsmittel := sample(c('Auto', 'Fahrrad'), size = .N, replace = TRUE)]

# Wähle Parteien basierend auf Einkommen, Beruf, Verkehrsmittel, Familienstand und Alter
wähler[, partei := sample(c("Grüne Partei", "Liberale Partei", "Soziale Partei", "Konservative Partei"),
                          size = 100000,
                          replace = T)]
wähler[beruf == 'Studierende/r' & verkehrsmittel == 'Fahrrad', partei := "Grüne Partei"]
wähler[beruf == 'Selbständig' & einkommen > 70000 & verkehrsmittel == 'Auto', partei := "Liberale Partei"]
wähler[beruf == 'Arbeiter/in' & verheiratet == 'Ja', partei := "Soziale Partei"]
wähler[alter > 50 & kinder > 1 & verkehrsmittel == "Auto", partei := "Konservative Partei"]

# Daten anzeigen
print(wähler)


# Stichprobe ziehen ----
stichprobe_n <- 100

umfrage <- wähler[sample(1:100000,
                         size = stichprobe_n,
                         replace = FALSE)]
print(umfrage)

umfrage[, table(partei)]

# plotten
ggplot(umfrage) +
    theme_minimal() +
    geom_histogram(aes(x = alter))

ggplot(umfrage) +
    theme_minimal() +
    geom_histogram(aes(x = partei), stat = "count")


# Viele Stichproben ziehen ----
set.seed(12345)

anzahl_umfrage <- 1000

umfrage_mehrfach <- rep_sample_n(wähler,
                                 size = stichprobe_n,
                                 replace = FALSE,
                                 reps = anzahl_umfrage)
umfrage_mehrfach = umfrage_mehrfach %>% data.table()

# große Variabilität
umfrage_anteile = umfrage_mehrfach[, .(anteil = .N / stichprobe_n), by = .(partei, replicate)][order(partei)]

# Stichprobenverteilungen
ggplot(data = umfrage_anteile,
       aes(x = anteil)) + 
    theme_minimal() +
    geom_histogram(binwidth = 0.01, boundary = 0.4, col = 'white') +
    facet_wrap(~ partei, scales = 'free_x') +
    labs(x = 'Anteile der Parteien', title = 'Verteilung der Parteien', y = 'Häufigkeit')


# Bootstrapping ----

## Stichprobe ziehen ----
stichprobe_n <- 1000

umfrage <- wähler[sample(1:100000,
                         size = stichprobe_n,
                         replace = FALSE)]

# wahrer Mittelwert des Alters in der Grundgesamtheit
wähler_mittelwert <- wähler[, mean(alter)]
wähler_mittelwert

## Bootstrap-Stichproben ziehen ----
stichprobe_n <- 1000
anzahl_umfrage <- 10000

bootstrap_umfrage <- rep_sample_n(wähler,
                                  size = stichprobe_n,
                                  replace = TRUE, # wichtig!
                                  reps = anzahl_umfrage)
bootstrap_umfrage = bootstrap_umfrage %>% data.table()

# mit replace kommen einzelne Wähler häufiger vor
bootstrap_umfrage[, .N, by = .(replicate, person_id)][order(-N)]

# Mittelwert aus Bootstrap berechnen
bootstrap_mittelwert <- bootstrap_umfrage[, .(Mittelwert = mean(alter)), by = replicate]
bootstrap_mittelwert[, mean(Mittelwert)]

# Standardfehler aus Quantilen berechnen
bootstrap_momente <- bootstrap_mittelwert[, .(mean = mean(Mittelwert),
                                              sd = sd(Mittelwert),
                                              lower = quantile(Mittelwert, 0.025),
                                              upper = quantile(Mittelwert, 0.975))]

# Plotten
ggplot() +
    theme_minimal() +
    geom_histogram(data = bootstrap_mittelwert,
                   aes(x = Mittelwert),
                   binwidth = 0.01) +
    geom_vline(xintercept = wähler_mittelwert, col = 'red') +
    geom_vline(aes(xintercept = mean), data = bootstrap_momente, col = 'blue') +
    geom_vline(aes(xintercept = lower), data = bootstrap_momente, col = 'blue', linetype = 'dashed') +
    geom_vline(aes(xintercept = upper), data = bootstrap_momente, col = 'blue', linetype = 'dashed') +
    labs(x = 'Mittelwert des Alters', title = 'Bootstrap-Mittelwert des Alters', y = 'Häufigkeit')


# Bootstrapped Anteil der Parteien
bootstrap_anteile = bootstrap_umfrage[, .(anteil = .N / stichprobe_n), by = .(partei, replicate)]
bootstrap_anteile[, mean(anteil), by = partei]

# Wahre Anteile der Parteien
wähler_anteile = wähler[, .(anteil = .N / 100000), by = partei]

# Plotten
ggplot(data = bootstrap_anteile,
       aes(x = anteil)) + 
    theme_minimal() +
    geom_histogram(binwidth = 0.01, col = 'white') +
    geom_vline(data = wähler_anteile, aes(xintercept = anteil), col = 'red') +
    facet_wrap(~ partei, scales = 'free_x') +
    labs(x = 'Anteile der Parteien', title = 'Bootstrapped Verteilung der Parteien', y = 'Häufigkeit')