# R/optim.ml.R In mhurdle: Multiple Hurdle Tobit Models

```inner_prod <- function(x, y = NULL){
x <- drop(x)
if (is.null(y)) y <- x else y <- drop(y)
sum(x * y)
}

outer_prod <- function(x, y = NULL){
x <- drop(x)
if (is.null(y)) y <- x else y <- drop(y)
outer(x, y)
}

cross_prod <- function(x, y = NULL){
x <- drop(x)
if (is.null(y)) y <- x else y <- drop(y)
drop(crossprod(x, y))
}

UP <- function(B, dx, dg)
B + (inner_prod(dx, dg) + inner_prod(dg, cross_prod(B, dg))) *
outer_prod(dx) / inner_prod(dx, dg) ^ 2 -
(outer_prod(cross_prod(B, dg), dx) +
t(outer_prod(cross_prod(B, dg), dx))) / inner_prod(dx, dg)

get_hessian <- function(x) if(! is.null(attr(x, "hessian"))) - attr(x, "hessian") else NULL
get_value <- function(x) - sum(as.numeric(x))

optim.ml <- function(f, start, g = NULL, H = NULL,
direction = c("max", "min"),
method = c("bfgs", "bhhh", "nr"),
value_tol = 1E-05,#sqrt(.Machine\$double.eps),
step_tol = 1E-05,#sqrt(.Machine\$double.eps),
iter_lim = 100, trace = 1L, ...){

init_time <- proc.time()
direction <- match.arg(direction)
method <- match.arg(method)
iter <- 0
compute_hessian <- (method == "nr")
x <- unname(start)
K <- length(x)
lnL <- f(x, gradient = TRUE, hessian = compute_hessian, ...)
if (method == "bfgs") B <- diag(K);solve(get_opg(lnL))#diag(K)
value <- get_value(lnL)
if (trace > 0L) cat(paste("Initial value of lnL:", round(value, 2), "\n"))
iter <- iter + 1
if (method == "bfgs") d_x <- - cross_prod(B, gradient)
if (method == "bhhh") d_x <- - solve(get_opg(lnL), gradient)
if (method == "nr") d_x <- - solve(get_hessian(lnL), gradient)
step <- 1
previous_lnL <- lnL
lnL <- f(x + d_x, hessian = compute_hessian, ...)
increasing <- get_value(lnL) < get_value(previous_lnL)
if (increasing){
while(get_value(lnL) < get_value(previous_lnL)){
step <- step * 2
previous_lnL <- lnL
lnL <- f(x + step * d_x, hessian = compute_hessian, ...)
}
step <- step / 2
lnL <- previous_lnL
}
else{
while(get_value(lnL) > get_value(previous_lnL)){
step <- step / 2
print(step)
previous_lnL <- lnL
lnL <- f(x + step * d_x, hessian = compute_hessian, ...)
}
}
x <- x + step * d_x

if (method == "bfgs") B <- UP(B, step * d_x, d_g)
print(get_value(previous_lnL))
print(get_value(lnL))
value_crit <- 1 - get_value(previous_lnL) / get_value(lnL)
if (is.na(value_crit)) value_crit <- 10
print(get_value(lnL))
print(value_crit)

if (value_crit < value_tol) code <- 2L
if (step < step_tol) code <- 3L
if (iter >= iter_lim) code <- 4L
if (trace > 0L) cat(paste("iteration:", iter, "lnL = ", round(get_value(lnL), 3),
"step: ", round(step, 4), "crit:", round(gradient_crit, 3), "\n"))
if (trace > 1L){
cat("--------------------------------------------\n")
print(round(resdet,3))
cat("--------------------------------------------\n")
}

est_stat <- structure(list(method = method,
value_crit = value_crit,
step_crit = step,
iter_crit = iter,
code = code,
time = proc.time() - init_time),
class = "est_stat")

result <- structure(list(lnL = get_value(lnL), coefficients = x, contributions = as.numeric(lnL),
opg = get_opg(lnL), hessian = get_hessian(lnL),
est_stat = est_stat),
class = "ml_max")
}
result
}

print.est_stat <- function(x, ...){

s <- round(x\$time[1],0)
h <- s %/% 3600
s <- s - 3600 * h
m <- s %/% 60
s <- s - 60 * m
cat(paste(x\$method, "method\n"))
tstr <- paste(h, "h:", m, "m:", s, "s", sep="")
cat(paste(x\$iter_crit,"iterations,",tstr,"\n"))
msg <- switch(x\$code,
"1" = "gradient close to zero",
"2" = "successive function values within tolerance limits",
"3" = "last step couldn't find higher value",
"4" = "iteration limit exceeded"
)
cat(paste(msg, "\n"))
}

print.ml_max <- function(x, ...) print(x\$est_stat)

summary.ml_max <- function(x){
se <- sqrt(diag(solve(x\$opg)))
round(cbind(x\$coefficients, se, student = x\$coefficients / se))
}
```

## Try the mhurdle package in your browser

Any scripts or data that you put into this service are public.

mhurdle documentation built on Dec. 11, 2021, 9:21 a.m.