#######################################################################################################
# MARSSaic function
# Adds to marssMLE object:
# elements in output arg
# samp.size, num.params
#######################################################################################################
MARSSaic = function(MLEobj, output=c("AIC","AICc"), Options=list(nboot=1000, return.logL.star=FALSE,
silent=FALSE)) {
# Options$nboot is the number of bootstrap replicates to do
# mssm.model is a specified model (needs the structure and parameter elements of the list)
# output tells what output to produce; this can be a vector if multiple items should be returned
# AIC, AICc, AICbp, AICbb, AICi, boot.params
# silent is a flag to indicate whether the progress bar should be printed
if(!is.list(Options)){
msg = " Options argument must be passed in as a list.\n"
cat("\nErrors were caught in MARSSaic \n", msg, sep="")
stop("Stopped in MARSSaic() due to problem(s) with arguments.\n", call.=FALSE)
}
## Set options if some not passed in
if(is.null(Options[["nboot"]])) Options$nboot = 1000
if(is.null(Options[["return.logL.star"]])) Options$return.logL.star = FALSE
if(is.null(Options[["silent"]])) Options$silent = FALSE
tmp = is.marssMLE(MLEobj)
if(!isTRUE(tmp)) {
if(!Options$silent) cat(tmp)
stop("Stopped in MARSSaic() due to problem(s) with the MLE object passed in.\n", call.=FALSE)
}
if(is.null(MLEobj[["logLik"]])){
msg = " No log likelihood. This function expects a model fitted via maximum-likelihood.\n"
cat("\n","Errors were caught in MARSSaic \n", msg, sep="")
stop("Stopped in MARSSaic() due to problem(s) with the MLE object passed in.\n", call.=FALSE)
}
return.list=list()
## Some renaming for readability
model = MLEobj$marss
#kf = MLEobj$kf
loglike = MLEobj$logLik
##### AIC and AICc calculations
if("AIC" %in% output | "AICc" %in% output){
K = 0
for (elem in c("Z", "A", "B", "U", "x0", "R", "Q","V0" )) {
pars = dim(model$free[[elem]])[2]
K = K + pars
}
MLEobj$AIC = -2*loglike + 2*K
samp.size = sum(!is.na(model$data))
MLEobj$AICc = ifelse(samp.size > (K+1),-2*loglike + 2*K*(samp.size/(samp.size-K-1)),"NA, number of data points less than K+1")
MLEobj$samp.size = samp.size
MLEobj$num.params = K
}
##### AICbb & AICbp
method = c("AICbb","AICbp")
for (m in method[which(method %in% output)]) {
bootstrap.method = switch(m,
AICbb = "innovations",
AICbp = "parametric")
logL.star=0
drawProgressBar = FALSE #If the time library is not installed, no prog bar
if(!Options$silent) { #then we can draw a progress bar
cat(paste(m, "calculation in progress...\n"))
prev = progressBar()
drawProgressBar = TRUE
}
boot.params = MARSSboot(MLEobj, nboot=Options$nboot, output="parameters", sim=bootstrap.method,
param.gen="MLE", silent=TRUE)$boot.params
if((bootstrap.method == "parametric") && ("boot.params" %in% output)) MLEobj$boot.params = boot.params
for(i in 1:Options$nboot) {
boot.model = MARSSvectorizeparam( MLEobj, parvec=boot.params[,i] ) #boot.model is a MLEobj
logL.star[i] = MARSSkf(boot.model, only.logLik=TRUE)$logLik
if(drawProgressBar) prev = progressBar(i/Options$nboot,prev)
}
MLEobj[[m]] = -4*(sum(logL.star))/Options$nboot + 2*loglike # -2*model$loglik + 2*(1/N)*(-2)*sum(boot$Yloglike-model$logLik)
if(Options$return.logL.star == TRUE) {
tmp = paste(m, ".logL.star", sep="")
MLEobj[[tmp]] = logL.star
}
}
return(MLEobj)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.