Nothing
vgFit <- function (x, freq = NULL, breaks = NULL, paramStart = NULL,
startMethod = "Nelder-Mead", startValues = "SL",
method = "Nelder-Mead", hessian = FALSE,
plots = FALSE, printOut = FALSE,
controlBFGS = list(maxit = 200),
controlNM = list(maxit = 1000),
maxitNLM = 1500, ...) {
xName <- paste(deparse(substitute(x), 500), collapse = "\n")
if (!is.null(freq)) {
if (length(freq) != length(x)) {
stop("vectors x and freq are not of the same length")
}
x <- rep(x, freq)
}
x <- as.numeric(na.omit(x))
startInfo <- vgFitStart(x = x, breaks = breaks, startValues = startValues,
paramStart = paramStart,
startMethodSL = startMethod,
startMethodMoM = startMethod, ...)
paramStart <- startInfo$paramStart
svName <- startInfo$svName
breaks <- startInfo$breaks
empDens <- startInfo$empDens
midpoints <- startInfo$midpoints
llfunc <- function(logParam) {
vgC <- logParam[1]
sigma <- exp(logParam[2])
theta <- logParam[3]
nu <- exp(logParam[4])
## added by Christine, email 23/09/2009
adjust <- 0.001
if (sigma < 0 | (abs(sigma - 0) < adjust)) {
sigma <- adjust
}
if (nu < 0 | (abs(nu - 0) < adjust)) {
nu <- adjust
}
## end of addition
return(-sum(log(dvg(x = x, param = c(vgC, sigma, theta, nu)))))
}
output <- numeric(7)
ind <- 1:4
if (method == "BFGS") {
opOut <- optim(paramStart, llfunc, NULL, method = "BFGS",
hessian = hessian, control = controlBFGS, ...)
}
if (method == "Nelder-Mead") {
opOut <- optim(paramStart, llfunc, NULL, method = "Nelder-Mead",
hessian = hessian, control = controlNM, ...)
}
if (method == "nlm") {
ind <- c(2, 1, 5, 4)
opOut <- nlm(llfunc, paramStart, hessian = hessian,
iterlim = maxitNLM, ...)
}
param <- as.numeric(opOut[[ind[1]]])[1:4]
param[2] <- exp(param[2])
param[4] <- exp(param[4])
names(param) <- c("vgC", "sigma", "theta", "nu")
maxLik <- -(as.numeric(opOut[[ind[2]]]))
conv <- as.numeric(opOut[[ind[4]]])
iter <- as.numeric(opOut[[ind[3]]])[1]
paramStart <- c(paramStart[1], exp(paramStart[2]),
paramStart[3], exp(paramStart[4]))
fitResults <- list(param = param, maxLik = maxLik,
hessian = if (hessian) opOut$hessian else NULL,
method = method, conv = conv, iter = iter, obs = x,
obsName = xName, paramStart = paramStart,
svName = svName, startValues = startValues,
breaks = breaks, midpoints = midpoints,
empDens = empDens)
class(fitResults) <- "vgFit"
if (printOut)
print(fitResults, ...)
if (plots) {
plot.vgFit(fitResults, ...)
}
return(fitResults)
}
print.vgFit <- function (x, digits = max(3, getOption("digits") - 3), ...)
{
if (!inherits(x,"vgFit")) {
stop("Object must belong to class vgFit")
}
cat("\nData: ", x$obsName, "\n")
cat("Parameter estimates:\n")
print.default(format(x$param, digits = digits), print.gap = 2,
quote = FALSE)
cat("Likelihood: ", x$maxLik, "\n")
cat("Method: ", x$method, "\n")
cat("Convergence code: ", x$conv, "\n")
cat("Iterations: ", x$iter, "\n")
invisible(x)
}
plot.vgFit <- function (x, which = 1:4,
plotTitles = paste(c("Histogram of ",
"Log-Histogram of ", "Q-Q Plot of ", "P-P Plot of "),
x$obsName, sep = ""),
ask = prod(par("mfcol")) < length(which) &&
dev.interactive(), ...)
{
if (!inherits(x,"vgFit")){
stop("Object must belong to class vgFit")
}
if (ask) {
op <- par(ask = TRUE)
on.exit(par(op))
}
par(mar = c(6, 4, 4, 2) + 0.1)
show <- rep(FALSE, 4)
show[which] <- TRUE
param <- x$param
breaks <- x$breaks
empDens <- x$empDens
mipoints <- x$midpoints
obs <- x$obs
obsName <- x$obsName
vgDens <- function(obs) dvg(obs, param = param,
tolerance = .Machine$double.eps ^ 0.5)
logvgDens <- function(obs) {
log(dvg(obs, param = param, tolerance = .Machine$double.eps ^ 0.5))
}
ymax <- 1.06 * max(vgDens(seq(min(breaks), max(breaks),
length = length(obs))),
empDens, na.rm = TRUE)
if (show[1]) {#histogram
hist.default(obs, breaks, right = FALSE, freq = FALSE,
ylim = c(0, ymax), main = plotTitles[1], ...)
curve(vgDens, add = TRUE, col = 2, ...)
title(sub = paste("param = (", round(param[1], 3), ",",
round(param[2], 3), ",", round(param[3], 3), ",",
round(param[4], 3), ")", sep = ""))
}
if (show[2]) {#log histogram
logHist(obs, breaks, include.lowest = TRUE, right = FALSE,
main = paste(plotTitles[2], obsName), ...)
curve(logvgDens, add = TRUE,
ylab = NULL, xlab =, col = 2, ...)
title(sub = paste("param = (", round(param[1], 3), ",",
round(param[2], 3), ",", round(param[3], 3), ",",
round(param[4], 3), ")", sep = ""))
}
if (show[3]) {#qq
qqvg(obs, param = param, main = paste(plotTitles[3], obsName), ...)
}
if (show[4]) {#pp
ppvg(obs, param = param, main = paste(plotTitles[4], obsName), ...)
}
invisible()
}
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.