Nothing
"disfitqua" <-
function(x, f, objfun=c("rmse", "mad"),
init.lmr=NULL, init.para=NULL, type=NA,
ptransf= function(t) return(t),
pretransf=function(t) return(t), verbose=FALSE, ... ) {
if(length(f) != length(x)) {
warning("Lengths of arguments 'f' and 'x' are not equal, returning NULL")
return(NULL)
}
if(is.na(type)) {
warning("Distribution type can not be NA, returning NULL")
return(NULL)
}
objfun <- match.arg(objfun)
n.para <- dist.list(type=type) # parameter count of the named distribution
if(! is.null(init.lmr)) {
if(length(init.lmr$L1) == 1) init.lmr <- lmorph(init.lmr)
if(length(init.lmr$lambdas) < n.para) {
warning("Length of the initial L-moments in argument 'init.lmr' is ",
"less than the number of parameters for the specified ",
"distribution type, returning NULL")
return(NULL)
}
}
# Having initial parameters will have priority over having initial L-moments
# because the optimization is based on the parameters changing not the L-moments
if(is.null(init.para)) {
# The premise of the "guessing" at the parameters is based on the quantiles
# but for many of intend applications (risk analysis) quantiles below the
# median are not usually available.
if(is.null(init.lmr)) {
if(verbose) {
message("Starting parameters not given, attempting to setup ",
"from the quantiles")
}
if(length(x) < n.para) {
warning("Length the data values in argument 'x' is less than the ",
"number of parameters for the specified distribution ",
"type, returning NULL")
return(NULL)
}
init.lmr <- rep(NA, n.para)
init.lmr[1] <- x[1] # treat the smallest quantile as the mean, in many
# expected applications, x[1] would be the median
# The mean() will work just fine without the 3rd item being present.
if(n.para > 1) init.lmr[2] <- mean(x[2:3], na.rm=TRUE) - x[1]
if(n.para > 2) init.lmr[3:n.para] <- 0
# The L-moments are now formed into the lmomco object and inherently
# checked for validity. The disfitqua() function assumes untrimmed
# versions if the L-moments are hacked from the quantiles.
try( init.lmr <- vec2lmom(c(init.lmr)) )
if(is.na(init.lmr$lambdas[1])) {
stop("FATAL: The initial L-moments (untrimmed version) deduced ",
"(assumed) from the quantiles are invalid through the ",
"vec2lmom() function, try providing alternative L-moments ",
"through the 'init.lmr' function or providing initial ",
"parameters through the 'init.para' function.")
return(NULL)
}
}
# The initial L-moments are used for the initial parameter estimates. The
# L-moments could either be hacked at using the quantiles themselves or
# from the staring L-moments.
init.para <- NULL
try( init.para <- lmom2par(init.lmr, type=type) )
if(is.null(init.para)) {
stop("FATAL: Initial parameters computing from either the deduced ",
"(assumed) L-moments from the quantiles or given L-moments by ",
"argument 'init.lmr' do not produce valid parameters for the ",
"distribution 'type' that has been given.")
return(NULL)
}
}
if(verbose) message(" Initial parameters: ",
paste(init.para, collapse=", "), "\n")
para <- init.para$para
fncontrol <- list(f=f, x=x, type=type)
OBJfun <- function(hat, obs) { NULL }
if(objfun == "rmse") {
OBJfun <- function(hat, obs) { sqrt( sum( (obs - hat)^2)/length(obs) ) }
} else if(objfun == "mad") {
OBJfun <- function(hat, obs) { sum(abs(obs - hat) )/length(obs) }
} else {
stop("should not be here in logic")
}
"fn" <- function(mypara) {
newpara <- NULL
try( newpara <- list(para=mypara, type=fncontrol$type) )
if(is.null(newpara)) {
warning("FATAL: new parameters could not be estimated, ",
"stepping the returning error to optim() to Inf")
return(Inf)
}
newXs <- par2qua(fncontrol$f, newpara)
#print(length(fncontrol$x))
if(length(newXs) != length(fncontrol$x)) {
warning("NEAR FATAL?: length of new quantiles does not equal length of ",
"the available quantiles, stepping returning error to optim() ",
"to Inf")
return(Inf)
}
err <- OBJfun(newXs, fncontrol$x)
return(err)
}
rt <- NULL
try( rt <- optim(ptransf(para), fn, ...) )
if(is.null(rt)) {
message("optim() call returned NULL, try changing the initial ",
"parameters or L-moments")
return(NA)
}
fit.para <- vec2par(pretransf(rt$par), type=type)
if(verbose) message(" Optimized parameters: ",
paste(fit.para, collapse=", "), "\n")
fit.para$init.para <- para
fit.para$source <- "disfitqua"
fit.para$disfitqua <- rt # preserve the operations of optim
return(fit.para)
}
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.