#### Bayesianische Schätzer für babyboom-Datensatz ####

#### Vorarbeiten ####

# Datensatz einlesen
baby <- read.table("Daten/babyboom.dat", header = TRUE)

# Zwischenzeiten berechnen
x <- diff(c(0,baby$GebZeit))

# ML-Schätzer
lambdaML <- 1/mean(x)


#### Aufgabe 1 c) ####

# Priori-Parameter (vgl. Angabe)
a <- 0.5
b <- 0.5

# Posteriori-Parameter
n <- length(x) # Anzahl der Daten
aPost <- a + n
bPost <- b + sum(x)

# Posteriori zeichnen
s <- seq(0, 0.1, 0.001)
plot(s, dgamma(s, shape = aPost, rate = bPost), type = "l",
     main = "Posteriori", xlab = expression(lambda), ylab = expression(paste("p(", lambda , "| x)")))
abline(v = aPost/bPost, col = 3) # Posteriori-Erwartungswert
abline(v = (aPost-1)/bPost, col = 4) # Posteriori-Modus
abline(v = lambdaML, col = 2, lty= 2) # ML-Schätzer
legend("topright", c("Posteriori-Erwartungswert", "Posteriori-Modus", "ML-Schätzer"), col = c(3,4,2), lty = c(1,1,2))


#### Aufgabe d) ####

# Ziehen aus der Posteriori
N <- 500
lambdaPost <- rgamma(N, shape = aPost, rate = bPost)
hist(lambdaPost, xlim = c(0, 0.1))
d <- density(lambdaPost) 
plot(d, main = "geschätzte Posteriori", xlim = c(0, 0.1))

# Empirische Schätzer für lambda
lambdaMean <- mean(lambdaPost) # Posteriori-Mittelwert
lambdaMean
abline(v = lambdaMean, col = 3)

lambdaMod <- d$x[which.max(d$y)] # Posteriori-Modus
lambdaMod
abline(v = lambdaMod, col = 4)

lambdaMed <- median(lambdaPost) # Posteriori-Median
lambdaMed
abline(v = lambdaMed, col = 5)

abline(v = lambdaML, col = 2, lty=  2) # ML-Schätzer (zum Vergleich)
legend("topright", c("Posteriori-Erwartungswert", "Posteriori-Modus", "Posteriori-Median", "ML-Schätzer"),
       col = c(3,4,5,2), lty = c(1,1,1,2))