Nothing
"qua2ci.cov" <-
function(x,f, type=NULL, nsim=1000, interval=c("confidence", "none"), level=0.90,
tol=1E-6, asnorm=FALSE, altlmoms=NULL, flip=NULL, dimless=TRUE,
usefastlcov=TRUE, nmom=5, getsimlmom=FALSE, verbose=FALSE, ...) {
interval <- match.arg(interval)
if(is.null(type)) {
message("must specify an lmomco distribution abbreviation, returning NA")
return(NA)
}
if(! any(dist.list() == type)) {
warning("invalid lmomco distribution abbreviation given, returning NA")
return(NA)
}
if(interval == "none" & length(f) > 1) {
warning("function is not vectorized in f for interval 'none', using ",
"only the first value")
f <- f[1]
}
if(! check.fs(f)) return(NA) # invalid probability
if(! check.fs(level) | level == 1) { # invalid confidence level
warning("level is not in [0,1), returning NA")
return(NA)
}
# Notice that if either of the special extensions of alternative L-moments
# and (or) the data have been flipped that dimensionless simulation will
# take place.
if(! is.null(altlmoms) & ! dimless) {
warning("Altervative L-moments are given, so setting 'dimless=TRUE'")
dimless <- TRUE
}
npar <- dist.list(type=type)
if(nmom < npar) {
warning("number of parameters distribution 'type' is smaller ",
"than 'nmom' resetting 'nmom to' ", npar)
nmom <- npar
}
# Note right off the bat that the data are actually flipped.
if(! is.null(flip)) x <- flip - x # flip the data if indicated
xlmoms <- lmoms(x, nmom=nmom) # L-moments of the 'original' or
# flipped data
if(is.null(altlmoms)) altlmoms <- xlmoms # set altlmoms if not already
if(verbose) {
message(" Lambdas (alternate or of x) "); print(altlmoms$lambdas)
}
if(! are.lmom.valid(altlmoms)) {
warning(" L-moments (either computed on 'x') or from 'altlmoms' ",
"are invalid, returning NA")
return(NA)
}
parent <- lmom2par(altlmoms, type=type)
if(verbose) {
message(" Parent parameters "); print(parent$para)
}
# Multi-variate Normal simulation of the L-moments
if(dimless) { # Removal of dimensions and scale standardization are requested
newx <- (x - xlmoms$lambdas[1]) / xlmoms$lambdas[2] # standardize
if(usefastlcov) { # Lmoments package has a fast algorithm
lmrcv <- Lmoments::Lmomcov(newx, rmax=nmom)
} else {
lmrcv <- lmoms.cov(newx, nmom=nmom)
}
sLM <- NULL
fed <- 0
if(nmom >= 2) fed <- c(fed, 1)
if(nmom >= 3) fed <- c(fed, altlmoms$ratios[3:nmom])
try(sLM <- MASS::mvrnorm(n=nsim, fed, lmrcv, tol=tol))
if(is.null(sLM)) {
message("x=")
print(x)
message("lmrcv=")
print(lmrcv)
warning("returning NULL (try increasing the tol-erance)")
return(NULL)
}
MU <- altlmoms$lambdas[1] # desired mean, which could be the original mean
L2 <- altlmoms$lambdas[2] # desired L2, which again could be the original
# based on afore logic involving: if(is.null(altlmoms)) altlmoms <- xlmoms
# notice that the LAMBDAS are converted, ratios will be RECOMPUTED later by
# logic involving the standard operation using first two lambdas and the
# trailing ratios in the vector-to-L-moments function vec2lmom().
# vec2lmom(c(slams[1:2], slams[3:nmom]/slams[2]))
for(i in 2:nmom) sLM[,i] <- sLM[,i]*L2 # rescale
sLM[,1] <- sLM[,1]*L2 + MU # rescale and shift
} else {
if(usefastlcov) { # Lmoments package has a fast algorithm
lmrcv <- Lmoments::Lmomcov(x, rmax=nmom)
} else {
lmrcv <- lmoms.cov(x, nmom=nmom)
}
sLM <- NULL
try(sLM <- MASS::mvrnorm(n=nsim, altlmoms$lambdas, lmrcv, tol=tol))
if(is.null(sLM)) {
message("x=")
print(x)
message("lmrcv=")
print(lmrcv)
warning("returning NULL (try increasing the tol-erance)")
return(NULL)
}
}
if(verbose) {
message(" Var-covar matrix"); print(lmrcv)
message(" Summary of simulated lambdas "); print(summary(sLM))
}
if(getsimlmom) return(sLM)
# Confidence limit computation
ci <- c((1-level)/2, 1-(1-level)/2) # division by 2 because two tailed
n <- length(f)
zz <- matrix(nrow=n, ncol=8)
for(i in 1:n) {
ff <- f[i]
quasf <- sapply(1:nsim, function(i) {
para <- NULL
slams <- sLM[i,] # extract L-moments, then cast to lmomco style
lmr <- vec2lmom(c(slams[1:2], slams[3:nmom]/slams[2]))
if(! are.lmom.valid(lmr)) return(NA) # silent warning
para <- lmom2par(lmr, type=type, ...)
if(is.null(para)) return(NA)
qf <- NULL
if(is.null(flip)) {
try(qf <- par2qua( ff, para))
ifelse(is.null(qf), return(NA), return(qf))
} else {
try(qf <- par2qua(1-ff, para))
ifelse(is.null(qf), return(NA), return(flip - qf))
}
})
if(interval == "none") return(quasf) # NAs will be returned by having this
# conditional ahead of the stripping and subsequent statistics thereof
if(any(is.na(quasf))) {
num.na <- length(quasf[is.na(quasf)])
warning("at least one of 'quasf' is NA, stripping ", num.na,
" NA values (entries) and continuing")
quasf <- quasf[! is.na(quasf)] # This should be quite unusual
# except for distributions having thin slices of valid L-moment
# domain or simulations near the boundary of the distributions,
# say above the Tau4 of the Generalized Logistic distribution.
}
if(length(quasf) < 3) {
warning("critical failure, no quantiles to then estimate the distribution for")
return(NULL)
}
lmr <- lmoms(quasf, nmom=3) # the L-moments of the simulated quantiles
mu <- lmr$lambdas[1]; lscale <- lmr$lambdas[2]
ifelse(asnorm, z <- qnorm(ci, mean=mu, sd=sd(quasf)), # assume Normality
z <- dat2bernqua(ci, quasf)) # nonparametric
fit <- ifelse(is.null(flip), par2qua( ff, parent),
flip - par2qua(1-ff, parent))
zz[i,] <- c(ff, z[1], fit, z[2],
median(quasf), mu, var(quasf), lscale)
}
# A quick review if nsim is large enough is that the qua_med/qua_mean
# columns should have values similar to the fit. The question of symmetry
# of the quantile estimates can be deduced by qua_med and qua_mean
# being nearly identical numbers.
colnames(zz) <- c("nonexceed", "lwr", "fit", "upr",
"qua_med", "qua_mean", "qua_var", "qua_lam2")
zz <- as.data.frame(zz)
return(zz)
}
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.