#' @export
summary.evmOpt <- function(object, nsim = 1000, alpha = .050, ...){
nc <- sapply(object$data$D, ncol)
if (all(nc == 1)){
env <- unclass(qqevm(object, nsim = nsim, alpha = alpha))
env <- list(data=sort(env$dat), envelope=env$sim, Q=env$p)
}
else {
x <- object
x$data$y <- object$residuals # Standard exponential if GPD model true
x$threshold <- 0
x$coefficients <- c(0, 0) # phi not sigma, so 0 not 1
env <- unclass(qqevm(x, nsim=nsim, alpha=alpha))
env <- list(data=sort(env$dat), envelope=env$sim, Q=env$p)
}
co <- cbind(object$coefficients, object$se, object$coefficients / object$se)
dimnames(co) <- list(names(coef(object)), c("Value", "SE", "t"))
res <- list(model = object, coefficients=co, envelope = env,
nsim = nsim, alpha = alpha)
oldClass(res) <- "summary.evmOpt"
res
}
#' @export
print.summary.evmOpt <- function(x, digits = 3 , ...){
oldX <- x
co <- coef(x)
env <- x$envelope
nsim <- x$nsim
alpha <- x$alpha
x <- x$model
cat("Call: ")
print(x$call, ... )
cat("\n")
print(x$family, verbose=TRUE, ...)
if (is.null(x$penalty) | x$penalty=="none"){
cat("\nModel fit by maximum likelihood.\n")
}
else {
cat("\nModel fit by penalized maximum likelihood.\n")
}
if (x$conv == 0) conv <- TRUE
else conv <- FALSE
cat( "\nConvergence:\t\t")
cat(conv)
if (x$rate < 1){
cat( "\nThreshold:\t\t")
cat(format(unname(x$threshold), digits=digits, ...))
cat( "\nRate of excess:\t\t")
cat(format(x$rate, digits=digits, ...))
}
cat("\n\n")
wh <- t(c(x$loglik, x$ploglik, AIC(x)))
colnames(wh) <- c("Log lik.", "Penalized log lik.", "AIC")
rownames(wh) <- ""
print(wh, print.gap=2, quote=FALSE, justify="left")
cat( "\n\nCoefficients:\n" )
print.default(format(co, digits=digits, ...), print.gap=2, quote=FALSE)
cat( "\n" )
cat(paste(nsim, " simulated data sets compared against observed data QQ-plot.\n", sep = "" ))
cat( "Quantile of the observed MSE: ", signif( env$Q, ... ), "\n")
out <- sum(env$data < env$envelope[1, ] | env$data > env$envelope[2, ])
level <- paste(round( 100 - alpha*100 , ...), "%", sep = "")
perc <- round(out / length( env$data ) * 100, digits=digits, ...)
perc <- paste(perc, "%", sep = "" )
cat( paste( out, " observations (", perc, ") outside the ", level, " simulated envelope.\n" , sep=""))
invisible(oldX)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.