Nothing
##====================================================================
## Author: Y. Deville
##
## The black-box maximization of the log-likelhood is inspired
## 'fitdistr' of Brian Ripley (package MASS)
##
## This is a re-factored code for "Renouv"
##
## o Make use of external function 'makeFuns' 'checkDist'
##
## o Allow the use of special distributions such as 'Lomax'
## and 'maxlo'.
##
## o Use concentrated log-likelihood for the historical case(s).
## The likelihood is concentrated with respect to 'lambda'
##
##
## TODO
##
## Doctorize possible problems with the 'maxlo' distribution
## if some of the historical data are not compatible with the value
## of the shape parameter.
##
##====================================================================
Renouv <- function(x,
threshold = NULL,
effDuration = NULL,
distname.y = "exponential",
MAX.data = NULL,
MAX.effDuration = NULL,
OTS.data = NULL,
OTS.effDuration = NULL,
OTS.threshold = NULL,
fixed.par.y = NULL,
start.par.y = NULL,
force.start.H = FALSE,
numDeriv = TRUE,
trans.y = NULL,
jitter.KS = TRUE,
pct.conf = c(95, 70),
rl.prob = NULL,
prob.max = 1.0-1e-4,
pred.period = NULL,
suspend.warnings = TRUE,
control = list(maxit = 300, fnscale = -1),
control.H = list(maxit = 300, fnscale = -1),
trace = 0,
plot = TRUE,
label = "",
...) {
## This could be passed as argument, but only TRUE is possible
## with historical data, and TRUE seems better.
info.observed <- TRUE
## for numerical differentiation
## eps <- sqrt(.Machine$double.eps) ## seems too small (for 2-nd order diff)
eps <- 1e-6
mc <- match.call()
OT <- TRUE
if (is(x, "Rendata")) {
if(trace) cat("processing the 'Rendata' object 'x'\n\n")
vn <- x$info$varName
x.OT <- x$OTdata[ , vn]
if (is.null(effDuration)) effDuration <- x$OTinfo$effDuration
if (is.null(threshold)) threshold <- x$OTinfo$threshold
} else {
if (length(threshold) == 0 || is.na(threshold) )
stop("a valid 'threshold' must be given")
if (length(x) > 0) {
if (length(effDuration) == 0 || is.na(effDuration) ||
(effDuration < 0))
stop("a valid 'effDuration' must be given")
x.OT <- x
} else {
if (!is.null(effDuration)) {
warning("'x' is NULL. No OT data provided, so 'effDuration' ",
"will be ignored")
}
effDuration <- 0
if (!(distname.y %in% c("exp", "exponential", "gpd", "GPD")) ) {
stop("'x' of length 0 is only allowed when 'distname.y' is ",
"\"exp\", \"exponential\", or \"gpd\"")
}
x.OT <- numeric(0)
OT <- FALSE
}
}
## transformation is only allowed with OT data
if ( !OT && !is.null(trans.y) ) {
stop("using 'trans.y' is only possible when OT data are provided")
}
## check and make MAXdata info
if (!missing(MAX.data)) {
MAX <- makeMAXdata(x = x, data = MAX.data, effDuration = MAX.effDuration)
} else {
MAX <- makeMAXdata(x = x, effDuration = MAX.effDuration)
}
## check and make OTSdata info
if (!missing(OTS.data)) {
OTS <- makeOTSdata(x = x, data = OTS.data, threshold = OTS.threshold,
effDuration = OTS.effDuration)
} else {
OTS <- makeOTSdata(x = x, threshold = OTS.threshold,
effDuration = OTS.effDuration)
}
if ( OTS$flag && (any(OTS$threshold < threshold)) )
stop("OTS thresholds must be >= threshold")
##=========================================================================
## build transformation functions, if necessary
##=========================================================================
start.par.y <- unlist(start.par.y)
fixed.par.y <- unlist(fixed.par.y)
myTransFuns <- transFuns(trans.y = trans.y,
distname.y = distname.y)
## shortcut names
transFlag <- myTransFuns$transFlag
transfun <- myTransFuns$transfun
invtransfun <- myTransFuns$invtransfun
##=========================================================================
## CODING RULES
## names ending with ".OT"used for objects related to OT data
##=========================================================================
nb.x <- length(x.OT)
if (nb.x > 0) {
ind <- (x.OT > threshold) & !is.na(x.OT)
if (!transFlag) {
y.OT <- as.double(x.OT[ind]) - threshold
} else {
threshold.trans <- myTransFuns$transfun(threshold)
dth <- threshold.trans - threshold
y.OT <- myTransFuns$transfun(as.double(x.OT[ind])) - threshold.trans
}
nb.OT <- length(y.OT)
if (!nb.OT) stop("no data above threshold")
} else {
y.OT <- numeric(0)
nb.OT <- 0
}
if (trace) cat("Number of obs > threshold", nb.OT, "\n")
##=========================================================================
## 'scale.OT' is a rounded quantity used to scale the parameters
## during optimisation
##
## 'y.OT' should be divided by 'scale.OT'. Instead, scale parameters
## of distributions should be divided by 'scale.y' (a vector), since we do
## not know in general which parameters are scale parameters.
##=========================================================================
if (nb.OT) {
mexpon.OT <- floor(log(quantile(y.OT, prob = 0.75), base = 10))
} else {
if (MAX$flag) {
mexpon.OT <- floor(log(mean(unlist(MAX$data) - threshold), base = 10))
} else if (OTS$flag) {
mexpon.OT <- floor(log(mean(unlist(OTS$data) - threshold), base = 10))
}
}
if (mexpon.OT >= 2) {
scale.OT <- 10^mexpon.OT
} else scale.OT <- 1.0
##=========================================================================
## check the distribution name and transfer essential objects
## to the current env
##=========================================================================
myDist <- checkDist(distname.y = distname.y, scale.OT = scale.OT)
funname.y <- myDist$funname.y
distname.y <- myDist$distname.y
if (myDist$special.y) {
parnames.y <- myDist$parnames.y
scale.y <- myDist$scale.y
} else {
parnames.y <- c(names(fixed.par.y), names(start.par.y))
scale.y <- rep(1, length.out = length(parnames.y))
}
parnames.all <- c("lambda", parnames.y)
parnb.y <- length(parnames.y)
##=========================================================================
## build probability functions and find the characteristics
## of the parameters. The two objects 'p.y' and 'fixed.y' are computed here
##=========================================================================
myFuns <- makeFuns(funname.y = funname.y,
parnames.y = parnames.y,
fixed.par.y = fixed.par.y,
trace = 0)
p.y <- myFuns$p.y
pf.y <- myFuns$pf.y
fixed.y <- myFuns$fixed.y
fixed.all <- c(FALSE, fixed.y)
names(fixed.all) <- parnames.all
ind.est.y <- (1:parnb.y)[!fixed.y]
## clean and complete 'myFuns', if necessary
funs <- list(trans = myTransFuns$transFlag,
transfun = myTransFuns$transfun,
invtransfun = myTransFuns$invtransfun,
dfun.y = myFuns$dfun.y,
pfun.y = myFuns$pfun.y,
qfun.y = myFuns$qfun.y,
logf.y = myFuns$logf.y,
q.y = myFuns$q.y,
F.y = myFuns$F.y)
##=========================================================================
## Prepare historical data if needed
##=========================================================================
hist.MAX <- MAX$flag
if (hist.MAX) {
z.MAX <- unlist(MAX$data)
if ( any(z.MAX <= threshold) ) {
stop("all historical 'MAX' data must exceed the threshold")
}
w.MAX <- MAX$effDuration
block.MAX <- MAX$block
r.MAX <- MAX$r
nblock.MAX <- length(r.MAX)
} else {
z.MAX <- numeric(0)
threshold.max <- numeric(0)
}
hist.OTS <- OTS$flag
if (hist.OTS) {
z.OTS <- unlist(OTS$data)
if ( any(z.OTS <= threshold) ) {
stop("all historical 'OTS' data must exceed the threshold")
}
w.OTS <- OTS$effDuration
threshold.OTS <- OTS$threshold
block.OTS <- OTS$block
r.OTS <- OTS$r
nblock.OTS <- length(r.OTS)
} else {
z.OTS <- numeric(0)
threshold.OTS <- numeric(0)
}
##=========================================================================
## compute the degree of freedom
## For historical blocks, each observation is considered as a
## df, and each empty OTS block must be so.
##=========================================================================
nobs <- length(y.OT)
if (hist.MAX) nobs <- nobs + length(z.MAX)
if (hist.OTS) nobs <- nobs + length(z.OTS) + sum(r.OTS == 0)
##=========================================================================
## Normal case: the ordinary sample OT exists
##=========================================================================
if (OT) {
## Event rate estimation
lambda.hat <- nb.OT / effDuration
est.N <- lambda.hat
cov.N <- lambda.hat / effDuration
names(est.N) <- "lambda"
names(cov.N) <- "lambda"
if (trace) {
cat("o Estimate of the process rate (evt/bloc length)", est.N, "\n")
}
##=====================================================================
## ML estimation ignoring historical data (if any). The main goal
## is here to compute
##
## 'est.y' estimates for th y part including fixed parms if any
## 'cov0.y' covariance matrix for the y parameters (filled with
## zeros for fixed params).
##
## For some distributions, a special routine may be used.
##
## XXX
##
## The MiwExp2 distribution should be handled as a special distribution
## in the future.
##
##=====================================================================
opt0 <- NULL
specials <- c("exponential", "weibull", "log-normal",
"gpd", "GPD", "gamma", "lomax", "maxlo")
if (distname.y %in% specials) {
resML <- fML(y = y.OT,
distname.y = distname.y,
parnames.y = parnames.y,
fixed.y = fixed.y,
fixed.par.y,
info.observed = TRUE)
est.y <- resML$estimate
cov0.y <- resML$cov
logLik0 <- resML$logLik
opt0 <- NULL
} else {
if (distname.y %in% c("MixExp2", "mixexp2")) {
warning("With the 'mixexp2' ML estimation may fail to converge")
distname.y <- "mixexp2"
funname.y <- "mixexp2"
parnames.y <- c("prob1", "rate1", "delta")
scale.y <- c(1, 1/scale.OT, 1/scale.OT)
## compute initial values
if ( is.null(start.par.y) ) {
ini <- ini.mixexp2(y.OT, plot = FALSE)$estimate
start.par.y <- c(ini["prob1"], ini["rate1"],
ini["rate2"]- ini["rate1"])
names(start.par.y) <- parnames.y
}
}
optim0 <- TRUE
## Arbitrary distribution: perform a general ML estimation
loglik0 <- function(parms) {
logL <- sum(funs$logf.y(parm = parms, x = y.OT))
}
## let's go...
if (trace) {
cat("o Optimization\n")
cat(" initial values\n")
print(start.par.y)
}
## Scale the parameters.
if (!is.null(control$parscale)) {
pn <- parnames.y[!fixed.y]
cp <- scale.y[!fixed.y]
cpn <- names(control$parscale)
if (!all(cpn %in% pn))
stop("'control$parscale' must have names in", pn)
cp[cpn] <- control$parscale[cpn]
control$parscale <- cp
} else {
control$parscale <- scale.y[!fixed.y]
}
control$ndeps <- rep(eps, length(control$parscale))
if (suspend.warnings) opt.old <- options(warn = -1)
if (trace) {
cat(" parscale used in 'control'\n"); print(control$parscale)
cat(" fnscale used in 'control'\n"); print(control$fnscale)
}
parLower.y <- myDist$parLower.y
parUpper.y <- myDist$parUpper.y
## if (is.null(parLower.y) && is.null(parUpper.y)) {
opt0 <- optim(par = start.par.y,
fn = loglik0,
method = "BFGS",
hessian = !numDeriv,
control = control)
if (suspend.warnings) opt.old <- options(opt.old)
if (opt0$convergence != 0) {
print(opt0)
stop("convergence not reached in optimisation")
}
logLik0 <- opt0$value
if (trace) {
cat(" optim ended normally.\n")
print(start.par.y)
}
## Now build est.y and cov0.y
est.y <- rep(NA, length(parnames.y))
names(est.y) <- parnames.y
## Fill the two parts separately
est.y[fixed.y] <- fixed.par.y
est.y[!fixed.y] <- opt0$par
## modified on 2010-01-25
cov0.y <- matrix(0, nrow = parnb.y, ncol = parnb.y)
colnames(cov0.y) <- rownames(cov0.y) <- parnames.y
## This was added in versions > 0.5-2 due to repeated
## problems with hessian
if (numDeriv) {
opt0Hessian <- numDeriv::hessian(func = loglik0, x = opt0$par)
} else {
opt0Hessian <- opt0$hessian
}
if (trace) {
cat("Hessian for ordinary logL\n")
print(opt0Hessian)
}
invHess <- try(solve(opt0Hessian))
if (!inherits(invHess, "try-error")) {
eig <- try(eigen(-invHess, symmetric = TRUE), silent = TRUE)
if (!inherits(eig, "try-error")) {
eig <- eig$values
if (any(eig <= 0)) {
warning("hessian not negative definite (1-st optim). ",
"Confidence limits will be misleading")
}
cov0.y[ind.est.y, ind.est.y] <- -invHess
} else {
warning("the eigenvalues of hessian could not be computed ",
"(1-st optim)")
}
} else {
warning("hessian could not be inverted (1-st optim)")
}
}
if (trace) {
cat("o Estimated values / covariance for the exceedances part\n")
print(est.y)
print(cov0.y)
}
param.N <- "lambda"
} else {
## No OT data usr parIni and drop the rate
if (hist.MAX) {
est.all <- parIni.MAX(MAX,
threshold = threshold,
distname.y = distname.y)
est.y <- est.all[-1L]
} else if (hist.OTS) {
est.all <- parIni.OTS(OTS,
threshold = threshold,
distname.y = distname.y)
est.y <- est.all[-1L]
} else {
stop("no data provided for estimation. Use 'RenouvNoEst'")
}
est.N <- NULL
cov.N <- NULL
cov0.y <- NULL
opt0 <- NULL
logLik0 <- NULL
}
##=========================================================================
## When historical data are present, the global log-likelihood must
## be maximised with "optim"
##=========================================================================
if ( hist.MAX || hist.OTS ) {
if (hist.MAX) { ## homogenise with other obs.
if (!transFlag) zMod.MAX <- z.MAX - threshold
else zMod.MAX <- transfun(z.MAX) - threshold.trans
zrMod.MAX <- tapply(zMod.MAX, block.MAX, min)
if (trace) {
cat("\n Take into account MAX historical data of time-length",
sum(w.MAX), "units\n")
}
}
if (hist.OTS) { ## homogenise with other obs.
if (!transFlag) {
zMod.OTS <- z.OTS - threshold
thresholdMod.OTS <- threshold.OTS - threshold
} else {
zMod.OTS <- transfun(z.OTS) - threshold.trans
thresholdMod.OTS <- transfun(threshold.OTS) - threshold.trans
}
if (trace) {
cat("\n Take into account OTS historical data on time-length",
sum(w.OTS), "units \n")
}
}
##======================================================================
## General Log-Likelihood: its formal args are obtained by
## concatenating "lambda" and the vector of parameters for the "y"
## part (exceedances).
##======================================================================
loglik <- function(parms) {
lambda <- parms[1]
lw <- lambda * effDuration
logL <- 0
if (OT) {
logL <- logL + dpois(nb.OT, lambda = lw, log = TRUE)
logL <- logL + sum(funs$logf.y(parm = parms[-1], x = y.OT))
}
if (hist.MAX) {
lw.MAX <- lambda * w.MAX
logL <- logL + sum(r.MAX * log(lw.MAX)) -
sum(lw.MAX * ( 1 - funs$F.y(x = zrMod.MAX, parm = parms[-1]) ))
logL <- logL + sum(funs$logf.y(parm = parms[-1], x = zMod.MAX))
}
if (hist.OTS) {
S.OTS <- (1.0 - funs$F.y(x = thresholdMod.OTS, parm = parms[-1]))
lw.OTS <- lambda * w.OTS
logL <- logL - sum(lw.OTS * S.OTS)
if (sum(r.OTS) > 0) {
logL <- logL + sum(r.OTS * log(lw.OTS)) +
sum(funs$logf.y(parm = parms[-1], x = zMod.OTS))
}
}
logL
}
##======================================================================
## Function needed to compute 'lambda' at optimum
##======================================================================
lambda.hat <- function(parms.y) {
## compute number of obs and discounted durations by blocks
w.TOT <- effDuration
r.TOT <- nb.OT
if (hist.MAX) {
S.MAX <- ( 1.0 - funs$F.y(x = zrMod.MAX, parm = parms.y) )
w.TOT <- w.TOT + sum(w.MAX * S.MAX)
r.TOT <- r.TOT + sum(r.MAX)
}
if (hist.OTS) {
S.OTS <- ( 1.0 - funs$F.y(x = thresholdMod.OTS, parm = parms.y) )
w.TOT <- w.TOT + sum(w.OTS * S.OTS)
r.TOT <- r.TOT + sum(r.OTS)
}
lambda.hat <- r.TOT / w.TOT
}
##=====================================================================
## Concentrated version (with respect to lambda. Caution:
## lambda is not in 'parms' here!!!
## ====================================================================
loglikC <- function(parms.y) {
## compute the numbers of obs and discounted durations by blocks
w.TOT <- effDuration
r.TOT <- nb.OT
if (hist.MAX) {
S.MAX <- ( 1.0 - funs$F.y(x = zrMod.MAX, parm = parms.y) )
w.TOT <- w.TOT + sum(w.MAX * S.MAX)
r.TOT <- r.TOT + sum(r.MAX)
}
if (hist.OTS) {
S.OTS <- ( 1.0 - funs$F.y(x = thresholdMod.OTS, parm = parms.y) )
w.TOT <- w.TOT + sum(w.OTS * S.OTS)
r.TOT <- r.TOT + sum(r.OTS)
}
lambda.hat <- r.TOT / w.TOT
## Ordinary OT part caution all blocks OT are grouped as one
logL <- r.TOT * log(lambda.hat)
logL <- logL + sum(funs$logf.y(parm = parms.y, x = y.OT))
if (hist.MAX) {
logL <- logL + sum(funs$logf.y(parm = parms.y, x = zMod.MAX))
}
if (hist.OTS) {
if (sum(r.OTS) > 0) {
logL <- logL + sum(funs$logf.y(parm = parms.y, x = zMod.OTS))
}
}
logL
}
## let's go...
if (trace) cat("o Optimisation\n")
if (suspend.warnings) opt.old <- options(warn = -1)
if (force.start.H) {
if (is.null(start.par.y))
stop(paste("You must provide initial values in 'start.par.y'",
"when 'force.start.H' is TRUE"))
par.ini.y <- start.par.y
## names(par.ini.y) <- parnames.y
} else {
## use the estimation without historical data. Note that
## the elts need to be named
par.ini.y <- est.y
names(par.ini.y) <- parnames.y
par.ini.y <- par.ini.y[!fixed.y]
if ( is.na(loglikC(par.ini.y)) || !is.finite(loglikC(par.ini.y)) ) {
## GPD "Weibullian" case for gpd. "doctorise" the
## identified problem The initial values may fail to
## be acceptable
if ( (distname.y %in% c("gpd", "GPD") && par.ini.y["shape"] < 0) ||
(distname.y == "maxlo") ) {
## if one of the historical levels is over the
## quantile with prob. 0.995 of the distribution,
## try to give a larger value to the scale
## parameter.
warning(paste("initial parameters for the GPD distribution",
"lead to some values outside of support. These",
"values are modified"))
upmax <- funs$q.y(parm = est.y, p = 0.999)
## known doctorizable situations
if ( (hist.MAX || hist.OTS) &&
any(c(z.MAX, z.OTS, threshold.OTS) > upmax + threshold) ) {
if (trace) {
mmax <- max(c(z.MAX, z.OTS, threshold.OTS))
cat("mmax = ", mmax, "\n")
if (distname.y == "gpd" || distname.y == "GPD") {
uu <- threshold - par.ini.y["scale"] / par.ini.y["shape"]
} else {
uu <- threshold - par.ini.y["scale"]
}
cat(" Upper limit of the support (estimation phase 1)", uu,"\n")
}
if (trace) {
cat(" Old par.ini.y[\"scale\"]", par.ini.y["scale"],"\n")
}
par.ini.y["scale"] <- par.ini.y["scale"] *
( max(c(z.MAX, z.OTS, threshold.OTS)) - threshold ) / upmax
if (trace) {
cat(" New par.ini.y[\"scale\"]", par.ini.y["scale"],"\n")
}
}
if (trace) {
cat(" Initial values 2nd stage (modified)\n"); print(par.ini.y)
}
} else {
stop(paste("loglik is NA or infinite at the initial values.\n",
"Give correct initial values in 'start.par.y'",
"and set 'force.start.H' to TRUE"))
}
}
}
if (trace) {
cat("o Initial values for 2nd stage optimization\n")
print(par.ini.y)
}
## Scale the parameters.
if (!is.null(control.H$parscale)) {
pn <- parnames.y[!fixed.y]
cp <- scale.y[!fixed.y]
cpn <- names(control.H$parscale)
if (!all(cpn %in% pn))
stop("'control$parscale' must have names in", pn)
cp[cpn] <- control.H$parscale[cpn]
control.H$parscale <- cp
} else {
control.H$parscale <- scale.y[!fixed.y]
}
control.H$ndeps <- rep(eps, length(control.H$parscale))
if (trace) {
cat(" parscale used in 'control.H'\n")
print(control.H$parscale)
cat(" ndeps used in 'control.H'\n")
print(control.H$ndeps)
cat(" fnscale used in 'control.H'\n")
print(control.H$fnscale)
}
opt <- optim(par = par.ini.y,
fn = loglikC,
hessian = !numDeriv,
control = control.H)
if (suspend.warnings) opt.old <- options(opt.old)
if (opt$convergence != 0) {
print(opt)
stop("convergence not reached in optimisation")
}
logLikFun <- loglik
logLik <- opt$value
## prepare Hessian evaluation
estimate1 <- rep(NA, parnb.y)
names(estimate1) <- parnames.y
estimate1[fixed.y] <- fixed.par.y
estimate1[!fixed.y] <- opt$par
lambda1 <- lambda.hat(opt$par)
estimate <- c(lambda = lambda1, estimate1)
parms1 <- c(lambda1, opt$par)
if (numDeriv) {
## require(numDeriv)
optHessian <- numDeriv::hessian(func = loglik, x = parms1)
} else {
##stop("you must use 'numDeriv = TRUE' since concentrated
##log-likelihood is used")
##
## XXX the hessian computation seems sometimes impossible
## with numDeriv.
##
opt1 <- optim(par = parms1,
fn = loglik,
hessian = TRUE,
control = list(maxit = 10, fnscale = -1))
optHessian <- opt1$hessian
}
if (trace) {
cat(" Fixed parameters in historical optimization and hessian\n")
print(fixed.all)
cat(" Maximising parameter\n")
print(opt$par)
cat(" Maximised log-likelihood\n")
print(opt$value)
cat(" Hessian\n")
print(optHessian)
}
## XXX Caution: only parameters which are not fixed.
cov.all <- matrix(0, nrow = parnb.y + 1, ncol = parnb.y + 1)
colnames(cov.all) <- rownames(cov.all) <- parnames.all
cov.all[!fixed.all, !fixed.all] <- NA
invHess <- try(solve(optHessian), silent = TRUE)
if (!inherits(invHess, "try-error")) {
eig <- try(eigen(-invHess, symmetric = TRUE), silent = TRUE)
if (!inherits(eig, "try-error")) {
eig <- eig$values
if (any(eig <= 0)) {
warning("hessian not negative definite (2-nd optim). ",
"Confidence limits will be misleading")
}
cov.all[!fixed.all, !fixed.all] <- -invHess
} else {
warning("the eigenvalues of hessian could not be computed (2-nd optim)")
}
} else {
warning("hessian could not be inverted (2-nd optim)")
}
res <- list(call = mc,
x.OT = x.OT,
y.OT = y.OT,
nb.OT = nb.OT,
effDuration = effDuration,
threshold = threshold,
distname.y = distname.y,
p.y = p.y,
parnames.y = parnames.y,
fixed.y = fixed.y,
trans.y = trans.y,
est.N = est.N,
cov.N = cov.N,
est.y = est.y,
cov.y = cov0.y,
corr.y = cov2corr(cov0.y),
estimate = estimate,
fixed = fixed.all,
df = p.y + 1L,
nobs = nobs,
p = p.y + 1L,
opt0 = opt0,
opt = opt,
logLik0 = logLik0,
logLik = logLik, ## constants for the historical parts.
logLikFun = logLikFun,
sigma = sqrt(diag(cov.all)),
cov = cov.all,
corr = cov2corr(cov.all),
history.MAX = MAX,
history.OTS = OTS,
transFlag = transFlag,
funs = funs)
## For later use...
est.y <- estimate[-1]
names(est.y) <- parnames.y
ind <- 1 + ((1L):(parnb.y))
cov.y <- cov.all[ind, ind, drop = FALSE]
} else {
##=============================================================
## No historical data, hence OT data. Just rearrange results
##=============================================================
estimate <- c(est.N, est.y)
cov.all <- matrix(0, nrow = parnb.y + 1L, ncol = parnb.y + 1L)
colnames(cov.all) <- rownames(cov.all) <- c("lambda", parnames.y)
cov.all[1, 1] <- cov.N
ind <- (1L):parnb.y
cov.all[1+ind, 1+ind] <- cov0.y[ind, ind, drop = FALSE]
cov.y <- cov0.y
logLikFun <- function(parms) {
logL <- sum(funs$logf.y(parm = parms, x = y.OT))
}
res <- list(call = mc,
x.OT = x.OT,
y.OT = y.OT,
nb.OT = nb.OT,
effDuration = effDuration,
threshold = threshold,
distname.y = distname.y,
p.y = p.y,
parnames.y = parnames.y,
fixed.y = fixed.y,
trans.y = trans.y,
est.N = est.N,
cov.N = cov.N,
est.y = est.y,
cov.y = cov0.y,
corr.y = cov2corr(cov0.y),
estimate = estimate,
fixed = fixed.all,
df = p.y + 1L,
nobs = nobs,
p = p.y + 1L,
opt = opt0,
logLik = logLik0,
logLikFun = logLikFun,
sigma = sqrt(diag(cov.all)),
cov = cov.all,
corr = cov2corr(cov.all),
history.MAX = MAX,
history.OTS = OTS,
funs = funs,
transFlag = transFlag)
}
##=========================================================================
## Compute a return level table
##=========================================================================
if (is.null(rl.prob)) {
rl.prob <- .rl.prob
rl.prob <- rl.prob[rl.prob <= prob.max]
} else {
if (any(is.na(rl.prob))) stop("'rl.prob' values can not be NA")
if ( any(rl.prob <= 0.0) || any(rl.prob >= 1.0) ) {
stop("'rl.prob' values must be > 0 and < 1")
}
rl.prob <- sort(rl.prob)
}
if (is.null(pred.period)) {
rr <- 2
pred.period <- (10^rr)*c(0.1, 0.2, 0.5, 1:10)
} else {
if (any(is.na(pred.period))) stop("'pred.period' values can not be NA")
pred.period <- sort(pred.period)
}
rl.period <- 1 / estimate[1] / (1 - rl.prob)
rl.sort <- sort(c(rl.period, pred.period), index.return = TRUE)
rl.period <- rl.sort$x
ind <- !duplicated(rl.period)
rl.period <- rl.period[ind]
ind.pred <- rl.period %in% pred.period
## Use pct.conf
ret.lev <- predict.Renouv(object = res,
newdata = rl.period,
## cov.rate = FALSE,
level = round(pct.conf) / 100,
trace = 1,
eps = 1e-6)
res$pct.conf <- pct.conf
rownames(ret.lev) <- NULL
res[["ret.lev"]] <- ret.lev
res[["pred"]] <- as.data.frame(ret.lev[ind.pred, , drop = FALSE])
res[["infer.method"]] <- attr(res[["pred"]], "infer.method")
##======================================================================
## perform a Kolmogorov-Smirnov test Note that a mix positional
## matching and name matching in the call!!! This is because
## 'F.y' has 'parm' as first arg
##======================================================================
if (nb.OT > 0) {
if (jitter.KS) {
KS <- ks.test(OTjitter(y.OT, threshold = 0.0),
funs$F.y, parm = estimate[-1])
} else KS <- ks.test(y.OT, funs$F.y, parm = estimate[-1])
res$KS.test <- KS
}
##======================================================================
## Perform Bartlett's (or Moran's) test for exponentiality and the
## LR test against a GPD alternative
##======================================================================
if (nb.OT > 0) {
if (distname.y == "exponential") {
res$expon.test <- gofExp.test(x = y.OT)
res$LRExp.test <- LRExp.test(x = y.OT)
}
}
class(res) <- "Renouv"
##======================================================================
## Return a model description for block maxima when possible
##======================================================================
if (distname.y %in% c("exponential", "exp")) {
RGEV <- Ren2gumbel(object = res)
cov.RGEV <- attr(RGEV, "vcov")
est.RGEV <- RGEV
## for cleaner shows
attr(est.RGEV, "threshold") <- NULL
attr(est.RGEV, "vcov") <- attr(est.RGEV, "jacobian") <- NULL
res$MAX <- list(distname = "gumbel",
blockDuration = 1,
estimate = est.RGEV,
sigma = sqrt(diag(cov.RGEV)),
cov = cov.RGEV)
}
if (distname.y %in% c("gpd", "lomax", "maxlo", "GPD")) {
RGEV <- Ren2gev(object = res)
cov.RGEV <- attr(RGEV, "vcov")
est.RGEV <- RGEV
## for cleaner shows
attr(est.RGEV, "threshold") <- NULL
attr(est.RGEV, "vcov") <- attr(est.RGEV, "jacobian") <- NULL
res$MAX <- list(distname = "gev",
blockDuration = 1,
estimate = est.RGEV,
sigma = sqrt(diag(cov.RGEV)),
cov = cov.RGEV)
}
if (plot) {
p <- try(plot(res, label = label, ...))
if (inherits(p, "try-error")) {
warning("an error occured in calling 'plot' frow 'Renouv'.",
" Try using 'plot' on the result")
}
}
return(res)
}
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.