#### Hilfsfunktionen für Laplace-Approximation im Beta-Binomial Modell ####

# Hilfsfunktion Logarithmierte Posteriori für Beta-Binomial Modell (nicht normalisiert)
#' @param pi Wert von pi, für den die logarithmierte Posteriori berechnet werden soll
#' @param n Anzahl der Versuche insgesamt, bekannter Parameter der Binomial-Verteilung
#' @param x Anzahl der positiven Ausgänge
#' @param a,b Parameter der Beta-Priori
#' 
#' @return Wert der nicht normalisierten, logarithmierten Posteriori an pi
H <- function(pi, n, x, a, b)
{
  return( dbinom(x, size = n, prob = pi, log = TRUE) + dbeta(pi, a, b, log = TRUE) )
}


# Laplace Approximation für Beta-Binomial Modell
#' @param H Funktion zur Berechnung der nicht normalisierten, logarithmierten Posteriori
#' @param g Transformationsfunktion des Parameters (z.B. Odds)
#' @param n Anzahl der Versuche (gesamt)
#' @param x Anzahl der positiven Ausgänge
#' @param a,b Parameter der Beta-Priori
#' @param eps obere und untere Schranke für die Optimierung. Default-Wert: 1e-12
#' 
#' @return Wert der Laplace-Approximation
laplaceApprox <- function(H, g, n, x, a, b, eps = 1e-12)
{
  # Logarithmierte Posteriori + log(g(pi))
  Q <- function(pi, n, x, a, b)
  {
    return( log(g(pi)) + H(pi, n, x, a, b) )
  }
  
  ### Optimierung
  # Startwert: 0.5
  # list(fnscale = -1): Maximierung statt Minimierung
  # Hessian = T: 2. Ableitung wird ausgegeben
  # lower: "fast 0", upper: "fast 1"
  resQ <- optim(0.5, Q, control = list(fnscale = -1), hessian = T, method = "L-BFGS-B", 
                lower = eps, upper = 1-eps, n = n, x = x, a = a, b = b)
  resH <- optim(0.5, H, control = list(fnscale = -1), hessian = T, method = "L-BFGS-B",
                lower = eps, upper = 1-eps, n = n, x = x, a = a, b = b)
  
  # Optimale Werte
  piQ <- resQ$par
  piH <- resH$par
  
  # 2. Ableitung der Funktionen an optimalen Werten
  zQ <- resQ$hessian
  zH <- resH$hessian
  
  # Laplace-Approximation
  return( sqrt(zH/zQ) * exp(Q(piQ, n, x, a, b) - H(piH, n, x, a, b)) )
}