hyperbReg <- function (X, y, breaks = NULL, Theta.start = NULL, start.values = c("fn"), method = c("BFGS"), trace = FALSE, controlbfgs = list(maxit = 100), controlnm = list(maxit = 1000), ...) { xlab <- deparse(substitute(y)) if (!is.null(dim(y))) { if (min(dim(y)) == 1) y <- as.vector(y) else stop("y must be a vector") } n <- length(y) if (missing(X)) { X <- as.matrix(rep(1, n)) Theta.names <- "mu" }else { if (is.null(colnames(X))){ X <- as.matrix(X) Theta.names <- outer(deparse(substitute(X)), as.character(1:ncol(X)), paste, sep = ".") }else Theta.names <- colnames(X) } Theta.names <- c("pi","zeta","delta",Theta.names) m <- ncol(X) if (is.null(Theta.start)) { qrX <- qr(X) beta <- qr.coef(qrX, y) resids <- qr.resid(qrX, y) if (is.null(breaks)) { breaks <- hist(resids, plot = FALSE, right = FALSE)$breaks } midpoints <- hist(resids, breaks, plot = FALSE, right = FALSE)$mids emp.dens <- hist(resids, breaks, plot = FALSE, right = FALSE)$density emp.dens <- ifelse(!is.finite(log(emp.dens)), NA, emp.dens) max.index <- order(emp.dens, na.last = FALSE)[length(emp.dens)] hyp.nu <- as.numeric(midpoints[max.index]) hyp.mu <- mean(resids) hyp.delta <- sd(resids) hyp.pi <- (hyp.nu - hyp.mu)/hyp.delta hyp.zeta <- 1 + hyp.pi^2 Theta.start <- c(hyp.pi, log(hyp.zeta), log(hyp.delta), beta) } else { if (length(Theta.start) != (3 + m)) stop("3+ncol(X) != length(Theta.start)") } output.bfgs <- NULL output.nm <- NULL #cat("Theta.start ",Theta.start,"\n") llfunc <- function(Theta) { K.nu <- besselK(exp(Theta[2]), nu = 1) location <- as.vector(X %*% as.matrix(Theta[-(1:3)])) likel <- 0 for(i in 1:length(y)) { likel <- likel - log(dhyperb(y[i], c(Theta[1:3], location[i]), K.nu = K.nu, log.pars=TRUE)); } } if (method == "BFGS") { optim.call <- optim(Theta.start, fn = llfunc, NULL, method = "BFGS", hessian = TRUE, control = controlbfgs) optim.par <- optim.call$par optim.val <- -optim.call$value optim.hessian <- optim.call$hessian } if (method == "Nelder-Mead") { optim.call <- optim(Theta.start, fn = llfunc, NULL, method = "Nelder-Mead", hessian = TRUE, control = controlnm) optim.par <- optim.call$par optim.val <- -optim.call$value optim.hessian <- optim.call$hessian } if (FALSE){ info <- optim.hessian if (all(is.finite(info))) { qr.info <- qr(info) info.ok <- (qr.info$rank == length(optim.par)) } else info.ok <- FALSE if (info.ok) { se2 <- diag(solve(qr.info)) se <- sqrt(ifelse(se >= 0, se2, NA)) } else se <- rep(NA, length(optim.par)) } Theta <- list(Pi=optim.par[1], Zeta=exp(optim.par[2]), Delta=exp(optim.par[3]), Beta=optim.par[-(1:3)]) list(call = match.call(), LogLikelihood = optim.val, Theta = Theta, optim = optim.call) }