#### Normalverteilung mit bekanntem Mittelwert und unbekannter Varianz ####

# Hilfsfunktion: Dichte der Invers-Gamma-Verteilung
#' @param s Wert, an dem die Dichte ausgewertet werden soll
#' @param a,b Parameter der Invers-Gamma-Verteilung
#' 
#' @return Wert der Invers-Gamma-Dichte an s
dInvGamma <- function(s, a, b)
{
  return( b^a / gamma(a) * s^(-a-1) * exp(-b/s) )
}

#### Prioris ####

# Konjugierte Priori
prioriK <- function(sigma2, a, b)
{
  return( dInvGamma(sigma2, a = a, b = b) )
}

# Jeffreys' Priori
prioriJ <- function(sigma2)
{
  return( 1/sigma2 )
}

# vorgegebene Werte
a <- 1
b <- 1e-3

# Plotten
s <- seq(0.01, 5, 0.01) # zum Zeichnen
plot(s, prioriK(s, a, b), type = "l", col = 3, ylim = c(0,5), 
     main = "Priori-Verteilungen", xlab = expression(sigma^2), ylab = expression(p(sigma^2)))
lines(s, prioriJ(s),  col = 4)
legend("topright", c("Konjugierte Priori", "Jeffreys' Priori"), col = c(3,4), lty = 1)

#### Posterioris ####

# Posteriori für konjugierte Priori
posterioriK <- function(sigma2, a, b, x, mu)
{
  n <- length(x)
  aPost <- a + n/2
  bPost <- b + 1/2*sum((x-mu)^2)
  
  return( dInvGamma(sigma2, a = aPost, b = bPost) )
}

# Posteriori für Jeffreys' Priori
posterioriJ <- function(sigma2, x, mu)
{
  n <- length(x)
  return( dInvGamma(sigma2, a = n/2, b = 1/2*sum((x-mu)^2) ) )
}

# Vorgegebene Werte
mu <- 1
x <- rnorm(10, mean = mu, sd = 1)  # d.h. wahres sigma^2 = 1

# Plotten
plot(s,  posterioriK(s, a, b, x, mu), type = "l", col = 3, ylim = c(0,2), 
     main = "Posteriori-Verteilungen", xlab = expression(sigma^2), ylab = expression(paste("p(", sigma^2 , "| x", mu, ")")))
lines(s,  posterioriJ(s, x, mu),  col = 4)
legend("topright", c("Konjugierte Priori", "Jeffreys' Priori"), col = c(3,4), lty = 1)


#### Auswirkung von n auf Posteriori ####

# Werte wie vorher
mu <- 1
a <- 1
b <- 1e-3
s <- seq(0.01, 5, 0.01)

# Teste verschiedene n
for(n in c(1,2,5,10,20, 50, 100, 200))
{
  x <- rnorm(n, mean = mu, sd = 1)  # d.h. wahres sigma^2 = 1
  
  # Posterioris
  plot(s, posterioriK(s, a, b, x , mu), type = "l", col = 3, ylim = c(0,2), 
       main = paste("Posteriori für n =", n), xlab = expression(sigma^2), ylab = expression(paste("p(", sigma^2 , "| x", mu, ")")))
  lines(s, posterioriJ(s, x, mu),  col = 4)
  legend("topright", c("Konjugierte Priori", "Jeffreys' Priori"), col = c(3,4), lty = 1)
  
  Sys.sleep(0.5)
}


#### Auswirkung von a,b auf Priori/Posteriori bei konjugierter Priori ####

# Werte wie vorher
mu <- 1
n <- 10
x <- rnorm(n, mean = mu, sd = 1)  # d.h. wahres sigma^2 = 1

# zum Zeichnen
s1 <- seq(0.01, 2, 0.01)
s2 <- seq(0.01, 10, 0.01)

### a variabel, b fest
aVals <- c(0.01, 0.1, 0.5, 1, 2, 5, 10)
b <- 1e-3

par(mfrow = c(1,2))
# Prioris
matplot(s1, sapply(aVals, prioriK, sigma2=s1, b=b), type = "l",  lty = 1, ylim = c(0,2), col = rainbow(length(aVals)),
        main = "Priori\nb = 1e-3", xlab = expression(sigma^2), ylab = expression(p(sigma^2)))
legend("topright", paste("a = ", aVals), col = rainbow(length(aVals)), lty = 1)

# Posterioris
matplot(s2, sapply(aVals, posterioriK, sigma2=s2, b=b, x=x , mu=mu),  type = "l", lty = 1, ylim = c(0,2),  col = rainbow(length(aVals)),
        main = "Posteriori\nb = 1e-3", xlab = expression(sigma^2), ylab = expression(paste("p(", sigma^2 , "| x", mu, ")")))
legend("topright", paste("a = ", aVals), col = rainbow(length(aVals)), lty = 1)


### a fest, b variabel
a <- 1
bVals <- c(1e-5, 1e-3, 1e-2, 0.1 , 1, 5, 10)

par(mfrow = c(1,2))
# Prioris
matplot(s1, sapply(bVals, prioriK, s=s1, a=a), type = "l",  lty = 1, ylim = c(0,2), col = rainbow(length(bVals)),
        main = "Priori\na = 1", xlab = expression(sigma^2), ylab = expression(p(sigma^2)))
legend("topright", paste("b = ", bVals), col = rainbow(length(bVals)), lty = 1)

# Posterioris
matplot(s2, sapply(bVals, posterioriK, s=s2, a=a, x=x , mu=mu),  type = "l", lty = 1, ylim = c(0,2),  col = rainbow(length(bVals)),
        main = "Posteriori\na = 1", xlab = expression(sigma^2), ylab = expression(paste("p(", sigma^2 , "| x", mu, ")")))
legend("topright", paste("b = ", bVals), col = rainbow(length(bVals)), lty = 1)