Skip to content

Commit

Permalink
Included support for the bgw package, which required fixed parameters…
Browse files Browse the repository at this point in the history 
… to be passed through the estimation environment as opposed to the functions.
  • Loading branch information
@edsandorf
edsandorf committed Jul 26, 2023
1 parent c0daea1 commit 6180bdf
Show file tree
Hide file tree
Showing 20 changed files with 247 additions and 161 deletions.
1 change: 1 addition & 0 deletions DESCRIPTION
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,7 @@ Imports:
maxLik,
numDeriv,
ucminf,
bgw,
generics
Suggests:
testthat,
Expand Down
1 change: 1 addition & 0 deletions NAMESPACE
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,6 +42,7 @@ export(get_param_free)
export(get_param_start)
export(get_v_chosen)
export(glance)
export(hessian_condition_number)
export(identify_choice_patterns)
export(inspect_list)
export(is_used)
Expand Down
80 changes: 57 additions & 23 deletions R/estimate-model.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,16 +30,13 @@ estimate_model <- function(type,
model <- tryCatch(switch(type,
maxlik = estimate_maxlik(log_lik,
param_free,
param_fixed,
workers,
ll,
control),
ucminf = estimate_ucminf(log_lik,
param_free,
param_fixed,
workers,
ll,
control)),
control),
bgw = estimate_bgw(log_lik,
param_free,
control)),
error = function(e) {
return(NULL)
})
Expand Down Expand Up @@ -76,18 +73,10 @@ estimate_model <- function(type,
#' @return A custom model object with the results fo the estimation.
estimate_maxlik <- function(log_lik,
param_free,
param_fixed,
workers,
ll,
control) {

model_obj <- maxLik::maxLik(log_lik,
start = param_free,
param_fixed = param_fixed,
workers = workers,
ll = ll,
return_sum = FALSE,
pb = NULL,
method = control[["method"]],
finalHessian = FALSE,
tol = control[["tol"]],
Expand All @@ -111,6 +100,7 @@ estimate_maxlik <- function(log_lik,

} else {
model[["convergence"]] <- FALSE

}

return(model)
Expand All @@ -125,26 +115,20 @@ estimate_maxlik <- function(log_lik,
#' @return A custom model object with the results of the estimation
estimate_ucminf <- function(log_lik,
param_free,
param_fixed,
workers,
ll,
control) {

model_obj <- ucminf::ucminf(par = param_free,
fn = log_lik,
hessian = 0,
param_fixed = param_fixed,
workers = workers,
ll = ll,
return_sum = TRUE,
pb = NULL,
control = list(
grtol = control[["gradtol"]],
xtol = control[["steptol"]],
stepmax = control[["stepmax"]],
maxeval = control[["iterlim"]],
grad = control[["grad"]],
gradstep = control[["gradstep"]]
gradstep = control[["gradstep"]],
trace = 5
))

# Added a minus to make the fit calculations correct
Expand All @@ -156,9 +140,59 @@ estimate_ucminf <- function(log_lik,

if (model_obj[["convergence"]] %in% c(1, 2, 3, 4)) {
model[["convergence"]] <- TRUE

} else {
model[["convergence"]] <- FALSE

}

return(model)
}

#' Estimate the model using the 'bgw' package
#'
#' The function is a convenient wrapper around \code{\link{bgw_mle}}.
#'
#' @inheritParams estimate_model
#'
#' @return A custom model object with the results of the estimation
estimate_bgw <- function(log_lik,
param_free,
control) {

model_obj <- bgw::bgw_mle(
calcR = log_lik,
betaStart = param_free,
calcJ = NULL, # Change if using non-numeric Jacobian
bgw_settings = list(
printLevel = control[["print_level"]]
)
)

# Added a minus to make the fit calculations correct
model <- list()
model[["optimum"]] <- -model_obj[["maximum"]]
model[["message"]] <- model_obj[["message"]]
model[["param_free"]] <- model_obj[["estimate"]]
model[["convergence_code"]] <- model_obj[["code"]]
model[["vcov"]] <- model_obj[["varcovBGW"]]
model[["param_stop"]] <- model_obj[["betaStop"]]
model[["gradient"]] <- model_obj[["gradient"]]
model[["hessian_condition_number_bgw"]] <- model_obj[["vcHessianConditionNumber"]]
model[["final_ll_bgw"]] <- model_obj[["finalLL"]]


# This is dangerous, but I don't know what return codes are successful conv.
model[["convergence"]] <- TRUE

# if (model_obj[["code"]] %in% c(1, 2, 3, 4)) {
# model[["convergence"]] <- TRUE
#
# } else {
# model[["convergence"]] <- FALSE
#
# }


return(model)
}
84 changes: 37 additions & 47 deletions R/estimate.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -80,14 +80,14 @@ estimate <- function(ll,
# ESTIMATE ----
cli::cli_h2("Estimating")
time_start_estimate <- Sys.time()

# Define the sets of parameters
# model_options <- validate(model_options, db_names_extra)

param_start <- unlist(model_options[["param"]])
param_free <- param_start[!(names(param_start) %in% model_options[["fixed"]])]
# This is just for the model object. The actual parameters are loaded as
# part of the estimation environment
param_fixed <- param_start[model_options[["fixed"]]]


model <- estimate_model(control[["optimizer"]],
log_lik,
param_start,
Expand All @@ -112,27 +112,18 @@ estimate <- function(ll,
sep = " "))

# OPTIMUM VALUES ----
model[["optimum_values"]] <- log_lik(get_param_free(model),
param_fixed,
workers,
ll,
return_sum = FALSE,
pb = NULL)
model[["optimum_values"]] <- log_lik(get_param_free(model))


# OPTIMUM AT 0 ----
ll_0 <- tryCatch({
ll_0_tmp <- log_lik(get_param_free(model) * 0,
param_fixed = param_fixed,
workers = workers,
ll =ll,
return_sum = TRUE,
pb = NULL)
ll_0_tmp <- log_lik(get_param_free(model) * 0)

if (control[["optimizer"]] %in% c("ucminf")) {
-ll_0_tmp

} else {
ll_0_tmp
}
switch(control[["optimizer"]],
ucminf = -ll_0_tmp,
maxlik = ll_0_tmp,
bgw = log(ll_0_tmp)
)

}, error = function(e) {
cli::cli_alert_danger("Failed to calculate 'll(0)'")
Expand All @@ -142,28 +133,31 @@ estimate <- function(ll,
model[["optimum_at_zero"]] <- ll_0

# GRADIENT ----
model[["gradient"]] <- numDeriv::grad(log_lik,
get_param_free(model),
param_fixed = param_fixed,
workers = workers,
ll = ll,
return_sum = TRUE,
pb = NULL)
# We only derive the gradient if we are not using 'bgw'
if (control[["optimizer"]] %in% c("maxlik", "ucminf")) {
model[["gradient"]] <- numDeriv::grad(log_lik,
get_param_free(model),
return_sum = TRUE)
}

# JACOBIAN ----
if (control[["calculate_jacobian"]]) {
cli::cli_h2("Calculating the Jacobian (gradient observations)")
time_start_jacobian <- Sys.time()

model[["gradient_obs"]] <- numDeriv::jacobian(log_lik,
get_param_free(model),
param_fixed = param_fixed,
workers = workers,
ll = ll,
return_sum = FALSE,
pb = NULL,
method = "simple")


if (control[["optimizer"]] == "bgw") {
model[["gradient_obs"]] <- numDeriv::jacobian(log_lik,
get_param_free(model),
method = "simple",
return_log = TRUE)

} else {
model[["gradient_obs"]] <- numDeriv::jacobian(log_lik,
get_param_free(model),
method = "simple")

}

time_diff <- Sys.time() - time_start_jacobian
cli::cli_alert_info(paste("Jacobian calculation took",
round(time_diff, 2),
Expand All @@ -179,15 +173,13 @@ estimate <- function(ll,
if (control[["calculate_hessian"]]) {
time_start_hessian <- Sys.time()
cli::cli_h2("Calculating the Hessian matrix")

# Calculate Hessian using the 'numderiv' package
model[["hessian"]] <- hessian("numderiv",
log_lik,
get_param_free(model),
param_fixed,
workers,
ll,
control[["method"]])
control[["method"]],
control[["optimizer"]])

# If calculation failed, try different method for calculating the Hessian
if (!hessian_complete(model[["hessian"]])) {
Expand All @@ -203,10 +195,8 @@ estimate <- function(ll,
model[["hessian"]] <- hessian("maxlik",
log_lik,
get_param_free(model),
param_fixed,
workers,
ll,
control[["method"]])
control[["method"]],
control[["optimizer"]])
}

# If the Hessian still cannot be calculated end and return the model object
Expand Down
21 changes: 11 additions & 10 deletions R/hessian.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,14 +9,13 @@
#' @inheritParams estimate_model
#' @param method A method for the approximation of the Hessian. Only necessary
#' for the maxlik package
#' @param optimizer A character string giving the optimizer used
#'
hessian <- function(type,
log_lik,
param_free,
param_fixed,
workers,
ll,
method) {
method,
optimizer) {

# Define n_par
n_par <- length(param_free)
Expand All @@ -29,20 +28,22 @@ hessian <- function(type,
width = 80
)

# Return log if we are using 'bgw'
return_log <- FALSE
if (optimizer == "bgw") {
return_log = TRUE
}

# Choose method for calculating Hessian.
hessian <- tryCatch(switch(type,
numderiv = numDeriv::hessian(func = log_lik,
x = param_free,
param_fixed = param_fixed,
workers = workers,
ll = ll,
return_log = return_log,
return_sum = TRUE,
pb = pb),
maxlik = maxLik::maxLik(logLik = log_lik,
start = param_free,
param_fixed = param_fixed,
workers = workers,
ll = ll,
return_log = return_log,
return_sum = TRUE,
pb = pb,
print.level = 0,
Expand Down
Loading

0 comments on commit 6180bdf

Please sign in to comment.