R/expert_zibinomial.R

In sparktseung/LRMoE: Logit-Weighted Reduced Mixture-of-Experts

# zibinomial Expert Function
# 17 Functions need to be implemented
######################################################################################
# Initialize the parameters of the zibinomial
######################################################################################
zibinomial.params_init <- function(y){
  # Estimate all the parameters that needed for further calculation
  n_init = max(y) + 2
  mu = mean(y)
  variance = var(y)
  p_init = ((variance + mu^2)/mu -1) / (n_init - 1)
  p_zero_init = 1 - mu/(n_init * p_init)
  return( list(p_zero = p_zero_init, n = n_init, p = p_init) )
}

zibinomial.exposurize <- function(params, exposure){
  # Calculate the exposures
  return( params )
}

zibinomial.set_params <- function(params){
  # Check the parameters are valid for distribution
  return( params )
}
######################################################################################
# Calculate the log likelihood and initialize the penalty function
######################################################################################
zibinomial.expert_ll_exact <- function(y, params){
  # If all the observations are exact, calculate its corresponding likelihood
  return( ifelse(y == 0,
                 log(params[["p_zero"]] + (1 - params[["p_zero"]])*binomial.pdf(params, x = y)),
                 log(1-params[["p_zero"]]) + binomial.logpdf(params, x = y)) )
}

zibinomial.expert_ll_not_exact <- function(tl, tu, yl, yu, params){
  return( faster_zi_result(tl, tu, yl, yu, params, "binomial") )
}

zibinomial.penalty <- function(params, penalty_params) {
  # Return the penalty applied on the parameters.
  # Keep in mind to set the default penalty parameters
  return(0)
  # stop("There are no approprate penalty parameters applied, please check `zibinomial.penalty()`")
}

zibinomial.default_penalty <- function() {
  return(binomial.default_penalty())
}

######################################################################################
# dzibinomial, pzibinomial, qzibinomial and rzibinomial implementations.
######################################################################################
zibinomial.simulation <- function(params, n) {
  return( (1 - rbinom(n, 1, params[["p_zero"]])) * rbinom(n, size = params[["n"]], prob = params[["p"]]) )
}

zibinomial.mean <- function(params) {
  # Calculate the mean based on the params
  return( (1 - params[["p_zero"]])*params[["n"]] * params[["p"]] )
}

zibinomial.variance <- function(params) {
  # Calculate the variance based on the params
  return( (1 - params[["p_zero"]]) * params[["n"]] * params[["p"]] * (1 - params[["p"]]) + params[["p_zero"]]*(1 - params[["p_zero"]]) * zibinomial.mean(params)^2 )
}

zibinomial.logpdf <- function(params, x) {
  # Return the log pdf based on the input x
  return( ifelse(is.infinite(x), 0, binomial.logpdf(params, x)) )
}

zibinomial.pdf <- function(params, x) {
  # Return the pdf based on the input x
  return( ifelse(is.infinite(x), -Inf, binomial.pdf(params, x)) )
}

zibinomial.logcdf <- function(params, q) {
  # return the log cdf based on the input x
  return( ifelse(is.infinite(q), 0, binomial.logcdf(params, q)) )
}

zibinomial.cdf <- function(params, q) {
  # return the cdf based on the input x
  return( ifelse(is.infinite(q), 1, binomial.cdf(params, q)) )
}

zibinomial.quantile <- function(params, p) {
  # return the percentage points based on the value of p
  return( ifelse(params[["p_zero"]] >= p, 0, binomial.quantile(params, p - params[["p_zero"]])) )
}

######################################################################################
# E Step, M Step and EM Optimization steps.
######################################################################################
# zibinomial._EStep <- function() {
#   # Perform the E step
#   NULL
# }
#
# zibinomial._MStep <- function() {
#   # Perform the M step
#   NULL
# }
#
# zibinomial.compute_EM <- function() {
#   # Perform the EM optimization
#   NULL
# }

zibinomial.EM_exact <- function(expert_old, ye, exposure, z_e_obs, penalty, pen_params) {
  # Perform the EM optimization with exact observations

  p_old = expert_old$get_params()$p_zero

  tmp_exp = ExpertFunction$new("binomial",
                               list(n = expert_old$get_params()$n,
                                    p = expert_old$get_params()$p),
                               pen_params)

  expert_ll_pos = tmp_exp$ll_exact(ye)

  z_zero_e_obs = z_e_obs * EM_E_z_zero_obs(ye, p_old, expert_ll_pos)
  z_pos_e_obs = z_e_obs - z_zero_e_obs
  p_new = EM_M_zero(z_zero_e_obs, z_pos_e_obs, 0.0, 0.0, 0.0)

  tmp_update = binomial.EM_exact(tmp_exp, ye, exposure, z_pos_e_obs,
                                         penalty, pen_params)

  return(list(p_zero = p_new,
              n = tmp_update$n,
              p = tmp_update$p))
}

zibinomial.EM_notexact <- function(expert_old,
                                           tl, yl, yu, tu,
                                           exposure,
                                           z_e_obs, z_e_lat, k_e,
                                           penalty, pen_params) {
  # Perform the EM optimization with exact observations

  p_old = expert_old$get_params()$p_zero

  if(p_old > 0.999999){
    return(list(p_zero = p_old,
                n = expert_old$n,
                p = expert_old$p))
  }

  tmp_exp = ExpertFunction$new("binomial",
                               list(n = expert_old$get_params()$n,
                                    p = expert_old$get_params()$p),
                               pen_params)

  expert_ll = rep(-Inf, length(yl))
  expert_tn_bar = rep(-Inf, length(yl))

  for(i in 1:length(yl)){
    expert_expo = tmp_exp$exposurize(exposure[i])
    result_set = expert_expo$ll_not_exact(tl[i], yl[i], yu[i], tu[i])
    expert_ll[i] = result_set[["expert_ll"]]
    expert_tn_bar[i] = result_set[["expert_tn_bar"]]
  }

  z_zero_e_obs = z_e_obs * EM_E_z_zero_obs(yl, p_old, expert_ll)
  z_pos_e_obs = z_e_obs - z_zero_e_obs
  z_zero_e_lat = z_e_lat * EM_E_z_zero_lat(tl, p_old, expert_tn_bar)
  z_pos_e_lat = z_e_lat - z_zero_e_lat
  p_new = EM_M_zero(z_zero_e_obs, z_pos_e_obs, z_zero_e_lat, z_pos_e_lat, k_e)

  tmp_update = binomial.EM_notexact(expert_old = tmp_exp,
                                            tl = tl, yl = yl, yu = yu, tu = tu,
                                            exposure = exposure,
                                            z_e_obs = z_pos_e_obs, z_e_lat = z_pos_e_lat,
                                            k_e = k_e,
                                            penalty = penalty, pen_params = pen_params)

  return(list(p_zero = p_new,
              n = tmp_update$n,
              p = tmp_update$p))
}
######################################################################################
# Register the zibinomial at zzz.R to the ExpertLibrary Object (Examples included)
######################################################################################