#### HPD-Intervall berechnen über Optimierung ####

# Hilfsfunktion zur Berechnung der Intervallänge bei vorgegebenem unteren Quantil
#' @param u unteres Quantil
#' @param g Kredibilitätsniveau (=1-alpha)
#' @param a,b Parameter der Gamma-Verteilung
#' 
#' @return Länge des Kredibilitätsintervalls mit u als unterer Grenze
gammaIntervalLength <- function(u, g, a, b)
{
  # Länge = Differenz der Quantile
  l <- qgamma(g+u, shape = a, rate = b) - qgamma(u, shape = a, rate = b)
  
  return(l)
}

# Berechnung des HPD-Intervalls über Optimierung
#' @param g Kredibilitätsniveau (=1-alpha)
#' @param a,b: Parameter der Gamma-Verteilung
#' 
#' @return Eine Liste mit drei Einträgen:
#'  interval: Ein Vektor mit der unteren und der oberen Grenze des Intervalls
#'  length: Die Länge des Intervalls
#'  quantile: Das zugehörige Quantil
gammaHPD <- function(g, a, b)
{
  # Unteren Grenze finden, die Intervallänge minimiert
  res <- optimize(f = gammaIntervalLength, interval = c(0, 1-g), g = g, a = a, b = b)
  
  u <- res$minimum
  
  return(list(interval = c(qgamma(u, shape = a, rate = b), qgamma(g+u, shape = a, rate = b)),
              length = res$objective,
              quantile = u))  
}