R Under development (unstable) (2023-08-28 r85029 ucrt) -- "Unsuffered Consequences" Copyright (C) 2023 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > library(Renext) Loading required package: evd > > ##======================================================================= > ## In this test, a sample of independent exceedances Over a Threshold > ## is drawn at random from an exponential distribution with mean > ## theta. The corresponding model is fitted. > ## > ## We the add the information that the (high) level 'X.U' > ## was never observed over aperiod of duration 'w.U'. We control the > ## theoretical relation between estimates > ## > ## theta.hat.new = theta.hat.old > ## - y.U * (lambda.hat.new / lambda.hat.old - 1) > ## > ## where 'old' is for OT only and 'new' is for OT AND Unobserved. > ## > ## > ## Note that the predicted return level curve will be far away > ## from the observed points, since the unobserved level has a high > ## leverage on the estimated parameters. Actually the unobserved > ## duaration is more than the double of that of the main sample. > ## > ##========================================================================= > > set.seed(123) > > n <- 20 > theta <- 100 > x <- 35 + theta*rexp(n) > > x.U <- 100 > w.U <- 100 + 100*runif(1) > > fitExp <- list() > > opar <- par(mfrow = c(2, 1), mar = c(5, 5, 3, 1)) > > fitExp[[1]] <- + Renouv(x = x, + effDuration = 50, + distname.y = "exponential", + threshold = 35, + conf.pct = c(70, 95), + prob.max = 0.99995, + pred.period = c(10, 100, 1000), + trace = 1, + main = "\"exponential\" without x.U") Number of obs > threshold 20 o Estimate of the process rate (evt/bloc length) 0.4 o Estimated values / covariance for the exceedances part rate 0.01232783 rate rate 7.598772e-06 Special inference for the exponential case without history Warning messages: 1: In predict.Renouv(object = res, newdata = rl.period, level = round(pct.conf)/100, : uncertainty on the rate not taken into account yet in the exponential with no history case 2: In regularize.values(x, y, ties, missing(ties)) : collapsing to unique 'x' values 3: In plot.window(...) : "conf.pct" is not a graphical parameter 4: In plot.xy(xy, type, ...) : "conf.pct" is not a graphical parameter 5: In axis(side = side, at = at, labels = labels, ...) : "conf.pct" is not a graphical parameter 6: In axis(side = side, at = at, labels = labels, ...) : "conf.pct" is not a graphical parameter 7: In box(...) : "conf.pct" is not a graphical parameter 8: In title(...) : "conf.pct" is not a graphical parameter > > fitExp[[2]] <- + Renouv(x = x, + effDuration = 50, + OTS.data = numeric(0), ## Unobserved + OTS.effDuration = w.U, + OTS.threshold = x.U, + distname.y = "exponential", + threshold = 35, + conf.pct = c(70, 95), + prob.max = 0.99995, + pred.period = c(10, 100, 1000), + control.H = list(fnscale = -1, trace = 3), + trace = 1, + main = "\"exponential\" plus unobserved lev. x.U ") Number of obs > threshold 20 o Estimate of the process rate (evt/bloc length) 0.4 o Estimated values / covariance for the exceedances part rate 0.01232783 rate rate 7.598772e-06 Take into account OTS historical data on time-length 147.7796 units o Optimisation o Initial values for 2nd stage optimization rate 0.01232783 parscale used in 'control.H' [1] 1 ndeps used in 'control.H' [1] 1e-06 fnscale used in 'control.H' [1] -1 Nelder-Mead direct search function minimizer function value for initial parameters = 143.129275 Scaled convergence tolerance is 2.13279e-06 Stepsize computed as 0.001233 BUILD 2 143.129275 142.325159 EXTENSION 4 142.325159 141.283770 REFLECTION 6 141.283770 140.849672 LO-REDUCTION 8 140.849672 140.820887 HI-REDUCTION 10 140.820887 140.820881 HI-REDUCTION 12 140.820881 140.817369 HI-REDUCTION 14 140.818236 140.817369 HI-REDUCTION 16 140.817581 140.817369 HI-REDUCTION 18 140.817420 140.817369 HI-REDUCTION 20 140.817381 140.817369 HI-REDUCTION 22 140.817371 140.817369 Exiting from Nelder Mead minimizer 24 function evaluations used Fixed parameters in historical optimization and hessian lambda rate FALSE FALSE Maximising parameter rate 0.01941634 Maximised log-likelihood [1] -140.8174 Hessian [,1] [,2] [1,] -421.6553 2719.076 [2,] 2719.0761 -91543.274 Warning messages: 1: In regularize.values(x, y, ties, missing(ties)) : collapsing to unique 'x' values 2: In plot.window(...) : "conf.pct" is not a graphical parameter 3: In plot.xy(xy, type, ...) : "conf.pct" is not a graphical parameter 4: In axis(side = side, at = at, labels = labels, ...) : "conf.pct" is not a graphical parameter 5: In axis(side = side, at = at, labels = labels, ...) : "conf.pct" is not a graphical parameter 6: In box(...) : "conf.pct" is not a graphical parameter 7: In title(...) : "conf.pct" is not a graphical parameter > > par(opar) > > estimates <- + array(NA, + dim = c(2, length(fitExp[[1]]$estimate)), + dimnames = list(c("sans", "US"), + names(fitExp[[1]]$estimate))) > > for (i in 1:2) estimates[i, ] <- fitExp[[i]]$estimate > > Test <- 1/estimates["US", "rate"] - + 1/estimates["sans", "rate"] - + (x.U - 35) * (estimates["US", "lambda"]/estimates["sans", "lambda"] -1) > > ## L'erruer doit être peite, disons ne pas dépasser 1% > RelError <- abs(Test)*estimates["sans", "rate"] > > stopifnot(RelError < 1.0) > > proc.time() user system elapsed 0.25 0.06 0.29