#### Normalverteilung mit Mittelwert und Präzision als unbekannte Parameter ####


#### Aufgabe b) ####

# Funktion zur Berechnung der gemeinsamen Posteriori von mu und kappa (nicht normalisiert) 
#' @param mu,kappa Werte der Parameter for Mittelwert und Präzision, für die die
#'   Posteriori ausgewertet werde soll
#' @param y Vektor der beobachteten Daten
#' 
#' @return Wert der nicht normalisierten Posteriori  
postJoint <- function(mu, kappa, y)
{
  n <- length(y)
  
  return(kappa^(n/2 - 1) * exp(- kappa/2 * sum( (y - mu)^2 ))) 
}

# Beobachtungen für y
y <- c(15, 20)

# mu und kappa - Werte
muVals <- seq(10, 25, 0.1)
kappaVals <- seq(0.01, 0.5, 0.01)

# gemeinsame Posteriori auf dem Gitter von mu und kappa-Werten auswerten
z <- matrix(NA, nrow = length(muVals), ncol = length(kappaVals)) # Ergebnis-Matrix

for(i in 1:length(muVals))
{
  for(j in 1:length(kappaVals))
    z[i,j] <- postJoint(mu = muVals[i], kappa = kappaVals[j], y)
}

# gemeinsame Posteriori plotten
contour(muVals, kappaVals, z, xlab = expression(mu), ylab = expression(kappa), main = "Gemeinsame Posteriori (nicht normalisiert)")

# Alternative Plot-Funktion
filled.contour(muVals, kappaVals, z, nlevels=20, color.palette=heat.colors,
               xlab = expression(mu), ylab = expression(kappa), main = "Gemeinsame Posteriori (nicht normalisiert)")


#### Aufgabe d) ####

source("dNonCentralT.R") # Hilfsfunktion laden

# Parameter aus der konkreten Stichprobe y berechnen
n <- length(y)
m <-  mean(y)
s2 <- sum( (y- mean(y))^2/(n^2-n))
nu <- n-1

# Plotten
plot(muVals, dNonCentralT(muVals, m = m, s = sqrt(s2) , nu), type = "l",
     xlab = expression(mu), ylab = expression(paste("Marginale Posteriori für ", mu)))