Nothing
fit.mle <- function(phi.init, loglik, theta.tx, extra.tx, Data, psi.names = NULL)
{
start.time <- Sys.time()
stopifnot(!is.null(Data$n))
n <- Data$n
# Combine the two lists: theta.tx(phi) and extra.tx(phi)
# into a vector. If the user provided psi.names, use those for the
# variable names.
psi.tx <- function(phi)
{
theta <- theta.tx(phi)
psi1 <- unlist(theta)
psi2 <- unlist(extra.tx(theta))
psi <- c(psi1, psi2)
if (!is.null(psi.names)) names(psi) <- psi.names
return(psi)
}
optim.method <- getOption("optim.method")
optim.control <- getOption("optim.control")
if (is.null(optim.control$fnscale)) { optim.control$fnscale <- -1 }
if (is.null(optim.control$maxit)) { optim.control$maxit <- 1e5 }
optim.res <- optim(par = phi.init, fn = loglik, method = optim.method,
control = optim.control, hessian = TRUE, Data = Data)
phi.hat <- optim.res$par
theta.hat <- theta.tx(phi.hat)
xi.hat <- extra.tx(theta.hat)
psi.hat <- psi.tx(phi.hat)
V.phi <- -solve(optim.res$hessian)
J.tx <- jacobian(psi.tx, phi.hat)
V.psi <- J.tx %*% V.phi %*% t(J.tx)
rownames(V.psi) <- colnames(V.psi) <- names(psi.hat)
loglik.hat <- optim.res$value
qq <- length(unlist(theta.hat))
aic <- -2 * loglik.hat + 2*qq
aicc <--2 * loglik.hat + 2*qq*n / (n-qq-1)
bic <- -2 * loglik.hat + qq*log(n)
# Note: follow SAS NLMIXED, treat Est/SE as t with df = n
df <- n
se <- sqrt(diag(V.psi))
t.val <- psi.hat / se
p.val <- 2 * (1 - pt(abs(t.val), df = df))
gr <- J.tx %*% grad(loglik, x = phi.hat, Data = Data)
estimates <- cbind(psi.hat, se, t.val, p.val, gr)
colnames(estimates) <- c("Estimate", "SE", "t-val", "P(|t|>t-val)", "Gradient")
elapsed.sec <- as.numeric(Sys.time() - start.time, 'secs')
res <- list(estimates = estimates, loglik = loglik.hat, aic = aic,
aicc = aicc, bic = bic, vcov = V.psi, V.phi = V.phi, optim.res = optim.res, df = df,
qq = qq, description = "<Default>", theta.hat = theta.hat,
xi.hat = xi.hat, phi.init = phi.init, elapsed.sec = elapsed.sec)
class(res) <- "mle.fit"
return(res)
}
#' @name S3 methods for mle.fit objects
confint.mle.fit <- function(object, parm = NULL, level = 0.95, ...)
{
if (!is.null(parm)) {
stop("Only supported value of parm is currently NULL")
}
dim.theta <- length(unlist(object$theta.hat))
dim.xi <- length(unlist(object$xi.hat))
na <- rep(NA, dim.theta + dim.xi)
DF <- data.frame(Estimate = na, SE = na, Lower = na, Upper = na)
rownames(DF) <- rownames(object$estimates)
w <- -qt((1-level)/2, df = object$df)
psi.hat <- object$estimates[,1]
se.psi.hat <- object$estimates[,2]
DF$Lower <- psi.hat - w * se.psi.hat
DF$Upper <- psi.hat + w * se.psi.hat
DF$Estimate <- psi.hat
DF$SE <- se.psi.hat
res <- list(ci = DF, df = object$df, t.quantile = w, fit = object, level = level)
class(res) <- "mle.fit.ci"
return(res)
}
#' Print summary of an \code{mle.fit.ci} object
#'
#' See generic \code{print} function for usage.
#' @param x An \code{mle.fit.ci} object from \code{confint.fit.mle}.
#' @param ... Not currently used.
#' @name S3 methods for mle.fit.ci objects
print.mle.fit.ci <- function(x, ...)
{
fit.out <- x$fit
dim.theta <- length(unlist(fit.out$theta.hat))
dim.xi <- length(unlist(fit.out$xi.hat))
printf("--- Parameter CIs (level %f) ---\n", x$level)
idx <- 1:dim.theta
print(x$ci[idx,])
if (dim.xi > 0)
{
printf("--- Additional CIs (level %f) ---\n", x$level)
idx <- 1:dim.xi + dim.theta
print(x$ci[idx,])
}
printf("--\n")
printf("Degrees of freedom: %d\n", fit.out$df)
printf("t-quantile: %f\n", x$t.quantile)
}
#' @name S3 methods for mle.fit objects
print.mle.fit <- function(x, ...)
{
printf("Fit for model:\n")
printf("%s\n", x$description)
dim.theta <- length(unlist(x$theta.hat))
dim.xi <- length(unlist(x$xi.hat))
DF <- as.data.frame(x$estimates)
DF[,1] <- round(DF[,1], 4)
DF[,2] <- round(DF[,2], 4)
DF[,3] <- round(DF[,3], 4)
DF[,4] <- my.numerical.format(DF[,4])
DF[,5] <- my.numerical.format(DF[,5])
printf("--- Parameter Estimates ---\n")
idx <- 1:dim.theta
print(DF[idx,])
if (dim.xi > 0)
{
printf("--- Additional Estimates ---\n")
idx <- 1:dim.xi + dim.theta
print(DF[idx,])
}
msg <- x$optim.res$message
printf("--\n")
printf("Elapsed Sec: %0.2f ", x$elapsed.sec)
printf("Degrees of freedom: %d\n", x$df)
printf("LogLik: %0.4f ", x$loglik)
printf("AIC: %0.4f ", x$aic)
printf("AICC: %0.4f ", x$aicc)
printf("BIC: %0.4f\n", x$bic)
printf("Converged status: %d ", x$optim.res$convergence)
printf("Message: %s\n", ifelse(is.null(msg), "<None>", msg))
}
#' @name S3 methods for mle.fit objects
summary.mle.fit <- function(object, ...)
{
object
}
#' @name S3 methods for mle.fit objects
coef.mle.fit <- function(object, ...)
{
object$estimates[,"Estimate"]
}
#' @name S3 methods for mle.fit objects
logLik.mle.fit <- function(object, ...)
{
object$loglik
}
#' @name S3 methods for mle.fit objects
AIC.mle.fit <- function(object, ..., k=2)
{
logLik(object) + k*object$qq
}
#' @name S3 methods for mle.fit objects
BIC.mle.fit <- function(object, ...)
{
object$bic
}
#' @name S3 methods for mle.fit objects
vcov.mle.fit <- function(object, ...)
{
object$vcov
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.