Package 'Renext'

Description

This package proposes fits and diagnostics for the so-called méthode du renouvellement, an alternative to other "Peaks Over Threshold" (POT) methods. The méthode du renouvellement generalises the classical POT by allowing the excesses over the threshold to follow a probability distribution which can differ from the Generalised Pareto Distribution (GPD). Weibull or gamma excesses are sometimes preferred to GPD excesses. The special case of exponential excesses (which falls in the three families: GPD, Weibull and gamma) has a special interest since it allows exact inference for the (scalar) parameter and for the quantiles form OT data (only).

The package allows the joint use of possibly three kinds of data or information. The first kind is classical excesses, or "OT data". It can be completed with two kinds of data resulting from a temporal aggregation, as is often the case for historical data. Both types are optional, and concern periods or blocks that must not overlap nor cross the OT period.

• MAX data correspond to the case where one knows the $r$ largest observations over each block. The number $r$ may vary across blocks. This kind of data is often called '$r$ largest', or "$r$ Largest Order Statistics" ($r$ LOS).

• OTS data (for OT Supplementary data) correspond to the case where one knows for each block $b$ all the observations that exceeded a threshold $u_b$ which is greater (usually much greater) than the main threshold $u$. The number $r_b$ of such observations can be zero, in which case we may say that $u_b$ is an unobserved level. A threshold $u_b$ is sometimes called a perception threshold.

Historical data are often available in hydrology (e.g. for river flood discharges, for sea-levels or sea surges) and can concern large periods such as past centuries. An unobserved level can typically be related to a material benchmark.

Maximum likelihood estimation is made possible in this context of heterogeneous data. Inference is based on the asymptotic normality of parameter vector estimate and on linearisation ("delta method") for quantiles or parameter functions.

The package allows the use of "marked-process observations" data (datetime of event and level) where an interevent analysis can be useful. It also allows the event dates to be unknown and replaced by a much broader block indication, e.g. a year number. The key point is then that the "effective duration" (total duration of observation periods) is known. Event counts for blocks can be used to check the assumption of Poisson-distributed events.

The package development was initiated, directed and financed by the french Institut de Radioprotection et de Sûreté Nucléaire (IRSN). The package is a non-academic tool designed for applied analysis on case studies and investigations or comparisons on classical probabilistic models.

Details

The DESCRIPTION file:

Index of help topics:

Brest                   Surge heights at Brest
Brest.years             Surge heights at Brest partial data
Brest.years.missing     Years with missing periods in 'Brest.year'
                        dataset
CV2                     Squared Coefficient of Variation
CV2.test                CV2 test of exponentiality
Dunkerque               Surge heights at Dunkerque
EM.mixexp               Expectation-Maximisation for a mixture of
                        exponential distributions
GPD                     Generalised Pareto Distribution
Garonne                 Flow of the french river La Garonne
Hpoints                 Plotting positions for exponential return
                        levels
Jackson                 Jackson's statistic
Jackson.test            Jackson's test of exponentiality
LRExp                   Likelihood Ratio statistic for exponential vs.
                        GPD
LRExp.test              Likelihood Ratio test of exponentiality vs. GPD
LRGumbel                Likelihood Ratio statistic for Gumbel vs. GEV
LRGumbel.test           Likelihood Ratio test for the Gumbel
                        distribution
Lomax                   Lomax distribution
Maxlo                   'maxlo' distribution
MixExp2                 Mixture of two exponential distributions
NBlevy                  Negative Binomial Levy process
OT2MAX                  Temporal aggregation of a Marked Process
OTjitter                Add a small amount of noise to a numeric vector
PPplot                  Diagnostic plots for Renouv objects
RLlegend                Legend management for return level plots
RLpar                   Graphical parameters for Return Level plots
RLplot                  Return level plot
Ren2gev                 Translate a vector of coefficients from a
                        Renewal-POT model with Pareto excesses into a
                        vector of GEV parameters
Ren2gumbel              Translate a vector of coefficients from a
                        Renewal-POT model with exponential excesses to
                        a vector of Gumbel parameters
Renext-package          Renewal Method for Extreme Values Extrapolation
Renouv                  Fit a 'Renouvellement' model
RenouvNoEst             Define a 'renouvellement' model without
                        estimation
SLTW                    Shifted Left Truncated Weibull (SLTW)
                        distribution
SandT                   Compute empirical survivals (S) and return
                        periods (T)
anova.Renouv            Compute an analysis of deviance table for two
                        nested Renouv objects
barplotRenouv           Barplot for Renouv "Over Threshold" counts
expplot                 Classical "exponential distribution" plot
fGEV.MAX                Fit a GEV distribution from block maxima or r
                        largest order statistics using an aggregated
                        Renewal POT process
fGPD                    Fit a two-parameters Generalised Pareto
                        Distribution from a sample
fgamma                  ML estimation of the Gamma distribution
flomax                  ML estimation of the Lomax distribution
fmaxlo                  ML estimation of a 'maxlo' distribution
fweibull                ML estimation of classical Weibull distribution
gev2Ren                 Translate a vector of GEV parameters into
                        renewal model
gof.date                Goodness-of-fit for the distribution of dates
gofExp.test             Goodness-of-fit test for exponential
                        distribution
gumbel2Ren              Translate a vector of Gumbel parameters into a
                        vector of parameters for a renewal model
ini.mixexp2             Simple estimation for the mixture of two
                        exponential distributions
interevt                Interevents (or interarrivals) from events
                        dates
logLik.Renouv           Log-likelihood of a "Renouv" object
mom.mixexp2             Moment estimation for the mixture of two
                        exponential distributions
mom2par                 Parameters from moments
pGreenwood1             Probability that the Greenwood's statistic is
                        smaller than one
parDeriv                Derivation of probability functions with
                        respect to the parameters
parIni.MAX              Initial estimation of GPD parameters for an
                        aggregated renewal model
plot.Rendata            Plot a Rendata object
plot.Renouv             Plot an object of class "Renouv"
predict.Renouv          Compute return levels and confidence limits for
                        a "Renouv" object
qStat                   Quantiles of a test statistic
rRendata                Simulate a random RenData object
readXML                 Read data using an XML index file
roundPred               Round quantiles in a pseudo-prediction table
skip2noskip             Fix non-skipped periods from skipped ones
spacings                Methods computing spacings between Largest
                        Order Statistics
summary.Rendata         Summary and print methods for "Rendata" objects
summary.Renouv          Summary and print methods for "Renouv" objects
translude               Make translucient colors
vcov.Renouv             Variance-covariance matrix of the estimates of
                        a "Renouv" object
weibplot                Classical Weibull distribution plot

This package contains a function Renouv to fit "renouvellement" models.

Author(s)

Yves Deville <[email protected]>, Lise Bardet <[email protected]>

Maintainer: Yann Richet <[email protected]>

References

• Miquel J. (1984) Guide pratique d'estimation des probabilités de crues, Eyrolles (coll. EDF DER).

• Coles S. (2001) Introduction to Statistical Modelling of Extremes Values, Springer.

• Embrechts P., Klüppelberg C. and Mikosch T. (1997) Modelling Extremal Events for Insurance and Finance. Springer.

See Also

The CRAN packages evd, ismev, extRemes, POT.

Examples

## 'Garonne' data set
summary(Garonne)
plot(Garonne)

## Weibull excesses
fG <- Renouv(x = Garonne,
             threshold = 3000,
             distname.y = "weibull",
             main = "Weibull fit for 'Garonne'")

coef(fG)
vcov(fG)
summary(fG)
logLik(fG)
## Re-plot if needed
plot(fG)

## Classical 'predict' method with usual formal args 
predict(fG, newdata = c(100, 150, 200), level = c(0.8, 0.9))

Description

Compute an analysis of deviance table for two nested Renouv objects

Usage

## S3 method for class 'Renouv'
anova(object, object1, trace = 1L, ...)

Arguments

Details

Of special interest is the case when the distribution of the excesses used in object is exponential while object1 uses a two-parameters alternative in the GPD family. We know then that the convergence to the asymptotic distribution is slow, and a numerical approximation of the true distribution of the test statistic is used when possible, i.e. when the objects do not use MAX or OTS data and the number of exceedances is between 8 and 500.

Value

An object of class "anova" inheriting from class "data.frame".

Note

The deviance of the models can not be interpreted: only the difference of the deviance is meaningful.

See Also

anova, LRExp.test.

Examples

## test using historical data
fit1Exp <- Renouv(Garonne,  distname.y = "exponential", plot = FALSE)
fit1GPD <- Renouv(Garonne, distname.y = "GPD", plot = FALSE)
anova(fit1Exp, fit1GPD)

## test without using historical data
x <- Garonne$OTdata$Flow
dur <- Garonne$OTinfo$effDuration

fit2Exp <- Renouv(x,  threshold = 2700,  effDuration = dur,
                  distname.y = "exponential", plot = FALSE)
fit2GPD <- Renouv(x, threshold = 2700, effDuration = dur,
                  distname.y = "GPD", plot = FALSE)
anova(fit2Exp, fit2GPD)

Description

Barplot for "Over Threshold" counts in time blocks (usually years)

Usage

barplotRenouv(data,
               blockname = colnames(data)[1],
               varname = colnames(data)[2],
               threshold = quantile(data[, varname], 0.2),
               na.block = NULL,
               plot = TRUE,
               main = NULL, xlab = NULL, ylab = NULL,
               mono = FALSE,
               prob.theo = 0.999,
               ...)

Arguments

Details

Blocks described in the na.block are omitted in the determination of counts. The object given in the na.block is coerced to character and the same is done for values of block before comparing them to the na.block values. If block variable is of class factor with levels representing years (e.g. 1980, 1981, etc.) missing blocks can be specified either as c("1980", "1981") or as numeric c(1980, 1981).

For the chi-square test, counts for neighbouring frequency classes are grouped in order to reach a minimum frequency of 5 in each group. E.g. if we expect respectively 1.0, 3.8 and 7.0 blocks with frequency 0, 1 and 2 for events, the three counts are grouped in one group with frequency 1.0+3.8+7.0=11.8. Note that this strategy of grouping is not unique and is likely to weaken the power of the test. Before grouping, the higher class theoretical probability is computed as the probability to obtain a count equal to or greater than the max value.

Value

A list with the following objects.

For both tests, the statistic follows a chi-square distribution under the null hypothesis . The list of results contains the statistic statistic, the number of degrees of freedom df and the $p$-value p.value.

Note

The two tests: (over-)dispersion and chi-square have one-sided (upper tail) $p$-value. In other words, we do not intend to reject when statistics take "abnormally small" values, but only when abnormally large values are met.

Author(s)

Yves Deville

References

See Yagouti A., Abi-Zeid I., Ouarda, T.B.M.J. and B. Bobée (2001), Revue de processus ponctuels et synthèse de tests statistiques pour le choix d'un type de processus Revue des Sciences de l'Eau, 1, pp. 323-361.

See Also

plot.Rendata

Examples

## na.block influence for Brest data
opar <- par(mfrow = c(2, 2))

bp1 <- barplotRenouv(data = Brest.years, threshold = 30,
         main = "missing periods ignored")
bp2 <- barplotRenouv(data = Brest.years, threshold = 30,
         na.block = 1992, main = "1992 missing")
bp3 <- barplotRenouv(data = Brest.years, threshold = 30,
         na.block = 1991:1993, main ="1991:1993 missing")
bp4 <- barplotRenouv(data = Brest.years, threshold = 30,
         na.block = Brest.years.missing, main = "all missing periods")

par(opar)

## threshold influence
opar <- par(mfrow = c(2,2))

thresh <- c(30, 35, 40, 50)

for (i in 1:length(thresh)) {
  bp  <- barplotRenouv(data = Brest.years, threshold = thresh[i],
                   na.block = Brest.years.missing,
                   main = paste("threshold =", thresh[i], "cm at Brest"))
}
par(opar)

Description

Surge heights near high tide at Brest tide gauge station (France), detailed version

Usage

Brest

Format

The format is: List of 5

• $info : List of 6

  • $name : chr "Brest"

  • $shortLab : chr "Surge Heights at Brest (France)"

  • $longLab : chr "Surge Heights near high tide, Brest (France)"

  • $varName : chr "Surge"

  • $varShortLab : chr "Surge"

  • $varUnit : chr "cm"

• $describe : chr "High tide sea surge over 30 cm at Brest (France)..."

• $OTinfo : List of 4

  • $start : chr POSIXct[1:1], format: "1846-01-01"

  • $end : chr POSIXct[1:1], format: "2009-01-01"

  • $effDuration : num 148

  • $threshold : num 30

• $OTdata : 'data.frame': 1289 obs. of 2 variables:

  • $date : POSIXct[1:1289], format: "1846-01-14" "1846-01-21" ...

  • $Surge : num [1:1289] 36 60 46 40 33 ...

• $OTmissing : 'data.frame': 43 obs. of 3 variables:

  • $start : POSIXct[1:43], format: "1846-01-01" "1847-01-01" ...

  • $end : POSIXct[1:43], format: "1846-01-04" "1847-01-21" ...

  • $comment : chr [1:43] "" "" "" "" ...

  - attr(*, "class")= chr "Rendata"

Details

Data are provided as a list.

• info gives general information about the data

• OTinfo gives general information about the Over the Threshold part of data. The effective duration (effDuration element) is the total duration for the periods with effective measurements.

• OTdata give OT measurements

• OTmissing gives start and end of the missing periods for OT measurements.

Data come from hourly sea levels measured and predicted by the french Service Hydrogéographique et Océanographique de la Marine (SHOM). Observed sea levels are available as REFMAR data at the url https://data.shom.fr/. Data were processed (declustered) by IRSN in order to provide a series of independent surge heights at high tide. Surge height at high tide is defined as the difference between the observed and the predicted maximal sea levels near high tide. A correction was applied to account for trend in the sea-level over the observation period.

The effective duration given in years is defined up to a small fraction of year due to leap years and leap seconds.

Source

https://data.shom.fr/

Examples

str(Brest)
Brest$OTinfo$start
plot(Brest)

Description

Surge heights at Brest (France)

Usage

Brest.years

Format

A data frame with 954 observations on the following 2 variables.

year

Year e.g; 1980

Surge

Surge heights above the threshold of 30 cm.

Details

These data are a simplified version of Brest. For each surge event only the year is retained as timestamp. Years with missing periods are available as a vector Brest.years.missing.

This dataset is useful for testing since similar data are sometimes met in the analyses.

Examples

names(Brest.years)

Description

Years with missing periods in the 'Brest.years' dataset

Usage

Brest.years.missing

Format

The format is: int [1:49] 1846 1847 1852 1857 1858 1859 1860 1861 1862 1863 ...

Details

Vector of years containing missing periods in the Brest.years dataset. This years should be ignored when computing yearly statistics such as event rates, since time records are lost.

Examples

print(Brest.years.missing)

Description

Squared Coefficient of Variation.

Usage

CV2(x)

Arguments

Details

Compute the squared Coefficient of Variation of one or several samples provided as a numeric vector or matrix.

Value

Numeric vector of the squared coefficients of variation.

Note

The squared coefficient of variation is the ratio $S^2/\bar{X}^2$ where $\bar{X}$ and $S^2$ are the sample mean and the sample variance. The variance is computed using the sample size $n$ as denominator, rather than the usual $n-1$.

Examples

n <- 30; nSamp <- 500
X <- matrix(rexp(n * nSamp), nrow= nSamp, ncol = n)
W <- CV2(X)
plot(density(W), main = "CV2 of exponential samples")

Description

Test of exponentiality based on the squared coefficient of variation.

Usage

CV2.test(x, method = c("num", "sim", "asymp"), nSamp = 15000)

Arguments

Details

The distribution of $\textrm{CV}^2$ is that of Greenwood's statistic up to normalising constants. It approximately normal with expectation $1$ and standard deviation $2/\sqrt{n}$ for a large sample size n. Yet the convergence to the normal is known to be very slow.

Value

A list of test results.

Note

This test is sometimes referred to as Wilk's exponentiality test or as WE1 test. It works quite well for a Lomax alternative (i.e. GPD with shape $\xi >0$), and hence can be compared to Jackson's test and the Likelihood-Ratio (LR) test of exponentiality. However, this test has lower power that of the two others while having a comparable computation cost due to the evaluation of the Greenwood's statistic distribution.

Author(s)

Yves Deville

References

S. Ascher (1990) "A Survey of Tests for Exponentiality" Commun. Statist. Theory Methods, 19(5), pp. 1811-1525.

See Also

The function CV2 that computes the statistic and LRExp.test or Jackson.test for functions implementing comparable tests or exponentiality with the same arguments.

Examples

n <- 30; nSamp <- 500
X <- matrix(rexp(n * nSamp), nrow = nSamp, ncol = n)
pVals <- apply(X, 1, function(x) CV2.test(x)$p.value)
plot(pVals)  ## should be uniform on (0, 1)

Description

Surge heights near high tide at Dunkerque tide gauge station (France)

Usage

Dunkerque

Format

The format is: List of 7

• $info : List of 6

  • $name : chr "Dunkerque"

  • $shortLab : chr "Surge Heights at Dunkerque (France)"

  • $longLab : chr "Surge Heights near high tide, Dunkerque (France)"

  • $varName : chr "Surge"

  • $varShortLab : chr "Surge"

  • $varUnit : chr "cm"

• $describe : chr "High tide sea surge over 30 cm at Dunkerque... "

• $OTinfo : List of 4

  • $start : POSIXct[1:1], format: "1956-01-01"

  • $end : POSIXct[1:1], format: "2009-01-01"

  • $effDuration: num 38.8

  • $threshold : num "30"

• $OTdata : 'data.frame': 740 obs. of 3 variables:

  • $date : POSIXct[1:740], format: "1956-11-27" "1956-12-03" ...

  • $Surge : num [1:740] 67.9 30.9 51.8 30.8 39.8 ...

  • $comment : chr [1:740] "" "" "" "" ...

• $OTmissing : 'data.frame': 83 obs. of 3 variables:

  • $start : POSIXct[1:83], format: "1956-01-01" "1956-08-08" ...

  • $end : POSIXct[1:83], format: "1956-06-07" "1956-11-03" ...

  • $comment: chr [1:83] "" "" "" "" ...

• $MAXinfo : 'data.frame' : 1 obs. of 3 variables:

  • $start : POSIXct[1:1], format: "1706-01-01"

  • $end : POSIXct[1:1], format: "1956-01-01"

  • $duration : num 250

• $MAXdata :'data.frame': 1 obs. of 4 variables:

  • $block : int 1

  • $date : POSIXct[1:1], format: "1953-02-01"

  • $Surge : num 213

  • $comment : chr "1"

  - attr(*, "class")= chr "Rendata"

Details

See Brest and Garonne datasets with the same list structure.

An 'historical' surge of 213 cm was observed on 1953-02-01 and is considered by experts as having a return period of 250 years.

Examples

Dunkerque$info
plot(Dunkerque)

Description

Experimental function for Expectation-Maximisation (EM) estimation

Usage

EM.mixexp(x, m = 2)

Arguments

Details

The EM algorithm is very simple for exponential mixtures (as well as for many other mixture models).

According to a general feature of EM, this iterative method leads to successive estimates with increasing likelihood but which may converge to a local maximum of the likelihood.

Value

List with

Note

The estimation is done for expectation (inverse rates) but the estimate vector in the result contains rates for compatibility reasons (e.g with exponential).

Author(s)

Yves Deville

See Also

mom.mixexp2 and ini.mixexp2 for "cheap" estimators when m = 2.

Examples

set.seed(1234)
x <- rmixexp2(n = 100, prob1 = 0.5, rate2 = 4)
EM.mixexp(x) -> res
res$estimate
matplot(res$Theta, type = "l", lwd = 2,
        xlab = "iteration", ylab = "theta",
        main = "exponential inverse rates")

Description

Plot a vector using "exponential distribution" scales

Usage

expplot(x,
        plot.pos = "exp",
        rate = NULL,
        labels = NULL,
        mono = TRUE,
        ...)

Arguments

Details

This plot shows $-\log[1-F(x)]$ against $x$ where $F(x)$ at point $i$ is taken as $i/(n+1)$ if plot.pos is "exp", or as the "median rank" approximation $(i-0.3)/(n+0.4)$ if plot.pos is "med".

If the data in x is a sample from an exponential distribution, the points should be roughly aligned. However the largest order statistics have high sampling dispersion.

Note

The log scale for y is emulated via the construction of suitable graduations. So be careful when adding graphical material (points, etc) to this graph with functions of the "add to plot" family (points, lines, ...).

The ML estimate of the rate parameter is the inverse of the sample mean.

Author(s)

Yves Deville

See Also

The weibplot function for a classical "Weibull" plot. The interevt is useful to compute interevents (or "interarrivals") that should follow an exponential distribution in the homogeneous Poisson process context.

Examples

x <- rexp(200)
 expplot(x, rate = 1/mean(x), labels = "fitted")

Description

Fast Maximum Likelihood estimation of the Gamma distribution.

Usage

fgamma(x, check.loglik = FALSE)

Arguments

Details

The likelihood is concentrated with respect to the scale parameter. The concentrated log-likelihood is a strictly concave function of the shape parameter which can easily maximised numerically.

Value

A list with the following elements

Note

The distribution is fitted by using the scale parameter rather than rate (inverse of scale).

Author(s)

Yves Deville

See Also

GammaDist in the stats package.

Examples

set.seed(9876)
alpha <- 0.06
beta <- rexp(1)
n <- 30
x <- rgamma(n, shape = alpha, scale = beta)
fit <- fgamma(x, check.loglik = TRUE)

## compare with MASS results
if (require(MASS)) {
   fit.MASS <- fitdistr(x, densfun = "gamma")
   rate <- 1 / fit$estimate["scale"]
   est <- c(fit$estimate, rate = rate)
   der <- rate * rate ## derivative of rate w.r.t scale
   sdest <- c(fit$sd, rate = der * fit$sd["scale"])
   tab <- rbind(sprintf(" %10.8f ", est),
                sprintf("(%10.8f)", sdest))
   colnames(tab) <- c("shape", "scale", "rate")
   rownames(tab) <- c("est", "sd")
   noquote(tab)
}

Description

Fit a GEV distribution from block maxima or $r$ largest order statistics using an aggregated Renewal POT process.

Usage

fGEV.MAX(MAX.data, MAX.effDuration,
         MAX = NULL,
         control = list(maxit = 300, fnscale = -1),
         scaleData = TRUE,
         cov = TRUE,
         info.observed = TRUE,
         trace = 0)

Arguments

Details

The data are assumed to provide maxima or $r$ largest statistics arising from an aggregated renewal POT model with unknown event rate $\lambda$ and unknown two-parameter Generalised Pareto Distribution for the excesses. A threshold $u$ is fixed below the given data and the three unknown parameters lambda, scale and shape of the POT model are found by maximising the likelihood. Then a vector of the three parameters for the GEV distribution is computed by transformation. The covariance matrix and standard deviations are computed as well using the jacobian matrix of the transformation.

The maximisation is for the log-likelihood with the rate lambda concentrated out, so it is a two-parameter optimisation.

Value

A list

Caution

Whatever be the data, the log-likelihood is infinite (hence maximal) for any vector of GEV parameters with shape $< -1$ and postive scale. Hence the log-likelihood should be maximised with a constraint on the shape, while an unconstrained optimisation is used here. In practice, for numerical reasons, the estimate usually remains inside the shape > -1 region. An estimation leading to shape $< -1$ must be considered as meaningless. An estimation with shape $< -0.5$ should be considered with care.

Note

This function could get more arguments in the future.

Author(s)

Yves Deville

References

The Renext Computing Details document.

Prescott P. and Walden A.T. (1980) Maximum Likelihood Estimation of the Parameters of the Generalized Extreme-Value Distribution. Biometrika 67(3), 723-724.

See Also

Renouv.

Examples

##====================================================================
## block maxima: simulated data and comparison with  the 'fgev'
## function from the 'evd' package
##====================================================================
set.seed(1234)
u <- 10
nBlocks <- 30
nSim <- 100   ## number of samples 
Par <- array(NA, dim = c(nSim, 3, 2),
             dimnames = list(NULL, c("loc", "scale", "shape"), c("MAX", "evd")))
LL <- array(NA, dim = c(nSim, 2),
            dimnames = list(NULL, c("MAX", "evd")))

for (i in 1:nSim) {
  rd <- rRendata(threshold = u,
                 effDuration = 1,
                 lambda = 12,
                 MAX.effDuration = rep(1, nBlocks),
                 MAX.r = rep(1, nBlocks),
                 distname.y = "exp", par.y = c(rate = 1 / 100))

  MAX <- Renext:::makeMAXdata(rd)
  fit.MAX <- fGEV.MAX(MAX = MAX)
  fit.evd <- fgev(x = unlist(MAX$data))
  Par[i, , "MAX"] <- fit.MAX$estimate
  Par[i, , "evd"] <- fit.evd$estimate
  LL[i, "MAX"] <- fit.MAX$loglik
  LL[i, "evd"] <- logLik(fit.evd)
}

##====================================================================
## r largest: use 'ismev::rlarg.fit' on the venice data set.
## NB 'venice' is taken from the 'evd' package here.
##====================================================================
## Not run:  
require(ismev);
fit1 <- ismev::rlarg.fit(venice)

## transform data: each row is a block
MAX.data <- as.list(as.data.frame(t(venice)))
## remove the NA imposed by the rectangular matrix format
MAX.data <- lapply(MAX.data, function(x) x[!is.na(x)])
MAX.effDuration <- rep(1, length(MAX.data))

fit2 <- fGEV.MAX(MAX.data = MAX.data,
                 MAX.effDuration = MAX.effDuration)

## estimates
est <- cbind(ismev = fit1$mle, RenextLab = fit2$estimate)
print(est)
# covariance
covs <- array(dim = c(2, 3, 3),
              dimnames = list(c("ismev", "RenextLab"),
                colnames(fit2$cov), colnames(fit2$cov)))
                
covs["ismev", , ] <- fit1$cov
covs["RenextLab", , ] <- fit2$cov
print(covs)

## End(Not run)

Description

Fit a two-parameters Generalised Pareto Distribution from a sample.

Usage

fGPD(x,
     info.observed = TRUE,
     shapeMin = -0.8,
     dCV = 1e-04,
     cov = TRUE,
     trace = 0)

Arguments

Details

This function mainly relies on the flomax and fmaxlo functions. When CV is larger than 1.0 + dCV, a Lomax distribution is fitted by Maximum Likelihood using a concentrated log-likelihood. When instead CV is smaller than 1.0 - dCV, a maxlo distribution is fitted. Finally, when CV -1.0 has absolute value <= dCV, an exponential distribution is fitted. In all cases, the result is translated into a parameter vector for the GPD.

Note that when CV is close to 1.0, fitting a Lomax or a maxlo distribution can lead to problems because the estimated values of the shape and scale parameter are large, and because the concentrated log-likelihood is a flat function of the scale parameter.

Value

A list

Note

It may happen that the estimated shape parameter is < -0.5, in which case the expected information matrix can not be computed, nor does the covariance matrix and standard deviations. In this case, the cov and sd objects contain NA values. This problem can arise also when the shape parameter is greater than but close to the value -0.5. Even when info.observed is TRUE, the information matrix, covariance and standard deviations are set to NA.

When the true (unknown) value is is < -0.5, the regularity conditions required in the ML approximated inference do not hold.

The default value of info.observed was set to TRUE from version 3.0-1 because standard deviations obtained with this choice are usually better.

Author(s)

Yves Deville

References

See the Renext Computing Details document.

See Also

The fmaxlo and flomax functions.

Examples

## Not run: 
set.seed(123456)
n <- 500
ns <- 1000
xi <- runif(ns, min = -0.5, max = 0.5)
X <- matrix(nrow = n, ncol = ns)

for (i in 1:length(xi)) {
  Xi <- rgpd(n, scale = 1, shape = xi[i])
  X[ , i] <- Xi
  res1 <- fGPD(Xi)
  res2 <- try(fpot(Xi, threshold = 0.0))
  if (inherits(res2, "try-error")) {
    cat(res2, "\n")
    break
  }
  logLik1 <- res1$loglik; logLik2 <- logLik(res2)
  if (abs(logLik1 - logLik2) > 0.001) {
    cat(sprintf("i = %d, xi = %7.4f\n", i, xi[i]))
    mat <- rbind(c(res1$estimate[1:2], logLik = logLik1),
                 c(res2$estimate[1:2], logLik = logLik2))
    rownames(mat) <- c("fGPD", "fpot")
    print(mat)
  }
}

## End(Not run)

Description

Fast Maximum Likelihood estimation of the Lomax distribution.

Usage

flomax(x,
       info.observed = TRUE,
       plot = FALSE,
       scaleData = TRUE,
       cov = TRUE)

Arguments

Details

The likelihood is concentrated with respect to the shape parameter. This function is increasing for small values of the scale parameter $\beta$. For large $\beta$, the derivative of the concentrated log-likelihood tends to zero, and its sign is that of $(1 - \textrm{CV}^2)$ where $\textrm{CV}$ is the coefficient of variation, computed using $n$ as denominator in the formula for the standard deviation.

The ML estimate does not exist when the sample has a coefficient of variation CV less than 1 and it may fail to be found when CV is greater than yet close to 1.

Value

A list with the following elements

Note

The estimates are biased for small or medium sized sample. The bias is positive for the shape parameter, thus the estimated shape tends to be larger than the true unknown value.

Fitting a Lomax distribution to an exponential sample might lead to a divergence since the exponential is the limit of a Lomax distribution with large shape and large scale with constant ratio shape/scale. Fitting this distribution to a sample having a coefficient of variation smaller than 1 is not allowed since it should lead to divergence of the estimation.

The default value of info.observed was set to TRUE from version 3.0-1 because standard deviations obtained with this choice are usually better.

Author(s)

Yves Deville

References

J. del Castillo and J. Daoudi (2009) "Estimation of the Generalized Pareto Distribution", Statist. Probab. Lett. 79(5), pp. 684-688.

D.E. Giles, H. Feng & R.T. Godwin (2013) "On the Bias of the Maximum Likelihood Estimator for the Two-Parameter Lomax Distribution" Comm. Statist. Theory Methods. Vol. 42, n. 11, pp. 1934-1950.

N. Johnson, S. Kotz and N. Balakrishnan Continuous Univariate Distributions vol. 1, Wiley 1994.

See Also

Lomax for the Lomax distribution.

Examples

## generate sample
set.seed(1234)
n <- 200
alpha <- 2 + rexp(1)
beta <- 1 + rexp(1)
x <- rlomax(n, scale = beta, shape = alpha)
res <- flomax(x, plot = TRUE)

## compare with a GPD with shape 'xi' and scale 'sigma'
xi <- 1 / alpha; sigma <- beta * xi  
res.evd <- evd::fpot(x, threshold = 0, model = "gpd")
xi.evd <- res.evd$estimate["shape"]
sigma.evd <- res.evd$estimate["scale"]
beta.evd <- sigma.evd / xi.evd 
alpha.evd <- 1 / xi.evd
cbind(Renext = res$estimate, evd = c(alpha = alpha.evd, beta = beta.evd))

Description

Fast Maximum Likelihood estimation of a 'maxlo' distribution.

Usage

fmaxlo(x,
       shapeMin = 1.25,
       info.observed = TRUE,
       plot = FALSE,
       scaleData = TRUE,
       cov = TRUE)

Arguments

Details

The 'maxlo' likelihood is concentrated with respect to the shape parameter, thus the function to be maximised has only one one scalar argument: the scale parameter $\beta$. For large scale $\beta$, the derivative of the concentrated log-likelihood tends to zero, and its sign is that of $(\textrm{CV}^2-1)$ where $\textrm{CV}$ is the coefficient of variation, computed using $n$ as denominator in the formula for the standard deviation.

The ML estimate does not exist when the sample has a coefficient of variation CV greater than 1.0 and it may fail to be found when CV is smaller than yet close to 1.0.

The expected information matrix can be obtained by noticing that when the r.v. $Y$ follows the 'maxlo' distribution with shape $\alpha$ and scale $\beta$ the r.v $V:= 1/(1-Y/\beta)$ follows a Pareto distribution with minimum 1 and and shape parameter $\alpha$. The information matrix involves the second order moment of $V$.

The default value of info.observed was set to TRUE from version 3.0-1 because standard deviations obtained with this choice are usually better.

Value

A list with the following elements

Note

The name of the distribution hence also that of the fitting function are still experimental and might be changed.

Author(s)

Yves Deville

See Also

Maxlo for the description of the distribution.

Examples

## generate sample
set.seed(1234)
n <- 200
alpha <- 2 + rexp(1)
beta <- 1 + rexp(1)
x <- rmaxlo(n, scale = beta, shape = alpha)
res <- fmaxlo(x, plot = TRUE)

## compare with a GPD with shape 'xi' and scale 'sigma'
xi <- -1 / alpha; sigma <- -beta * xi
res.evd <- evd::fpot(x, threshold = 0, model = "gpd")
xi.evd <- res.evd$estimate["shape"]
sigma.evd <- res.evd$estimate["scale"]
beta.evd <- -sigma.evd / xi.evd 
alpha.evd <- -1 / xi.evd
cbind(Renext = res$estimate, evd = c(alpha = alpha.evd, beta = beta.evd))

Description

Fast Maximum Likelihood estimation of the classical two parameters Weibull distribution.

Usage

fweibull(x, info.observed = TRUE, scaleData = TRUE, cov = TRUE,
         check.loglik = FALSE)

Arguments

Details

The ML estimates are obtained thanks to a reparameterisation with $\eta = scale^{1/shape}$ in place of shape. This allows the maximisation of a one-dimensional likelihood $L$ since the $\eta$ parameter can be concentrated out of $L$. This also allows the determination of the expected information matrix for $[shape,\,\eta]$ rather than the usual observed information.

Value

A list

Note

The default value of info.observed was set to TRUE from version 3.0-1 because standard deviations obtained with this choice are usually better.

Author(s)

Yves Deville

See Also

weibplot for Weibull plots.

Examples

n <- 1000
set.seed(1234)
shape <- 2 * runif(1)
x <- 100 * rweibull(n, shape = 0.8, scale = 1)
res <- fweibull(x)

## compare with MASS
if (require(MASS)) {
   res2 <- fitdistr(x , "weibull")
   est <- cbind(res$estimate, res2$estimate)
   colnames(est) <- c("Renext", "MASS")
   loglik <- c(res$loglik, res2$loglik)
   est <- rbind(est, loglik)
   est
}

## Weibull plot
weibplot(x,
         shape = c(res$estimate["shape"], res2$estimate["shape"]),
         scale = c(res$estimate["scale"], res2$estimate["scale"]),
         labels = c("Renext 'fweibull'", "MASS 'fitdistr'"),
         mono = TRUE)

Description

Flow of the french river La Garonne at le Mas d'Agenais

Usage

Garonne

Format

The format is: List of 7

• $info : List of 6

  • $name : chr "Garonne"

  • $shortLab : chr "La Garonne at Le Mas d'Agenais"

  • $longLab : chr "River flow of La Garonne at Le Mas d'Agenais"

  • $varName : chr "Flow"

  • $varShortLab : chr "Flow"

  • $varUnit : chr "m3/s"

• $describe : chr "Flow of the french river La Garonne ..."

• $OTinfo :List of 4

  • $start : POSIXct[1:1], format: "1913-01-01"

  • $end : POSIXct[1:1], format: "1978-01-01"

  • $effduration : num 65

  • $threshold : num 2500

• $OTdata : 'data.frame': 151 obs. of 3 variables:

  • $date : POSIXct[1:151], format: "1913-04-08" "1913-04-25" ...

  • $Flow : num [1:151] 2600 2800 2700 4579 3400 ...

  • comment : chr [1:151] "" "" "" "" ...

• $OTmissing : NULL

• $MAXinfo :'data.frame': 1 obs. of 3 variables:

  • $start : POSIXct[1:1], format: "1770-01-01"

  • $end : POSIXct[1:1], format: "1913-01-01"

  • $duration : num 143

• $MAXdata :'data.frame': 12 obs. of 4 variables:

  • $block : num [1:12] 1 1 1 1 1 1 1 1 1 1 ...

  • date : POSIXct[1:12], format: NA NA ...

  • $Flow : num [1:12] 7500 7400 7000 7000 7000 6600 6500 6500 6400 6300 ...

  • $comment : chr [1:12] "1 (1875)" "2 (1770)" "3 (1783)" "4 (1855)" ...

  - attr(*, "class")= chr "Rendata"

Details

The data concern the french river La Garonne at the gauging station named Le Mas d'Agenais where many floods occurred during the past centuries.

The data consist in OT data and historical data. The variable is the river flow in cube meter per second $(\textrm{m}^3/\textrm{s})$ as estimated from the river level using a rating curve. The precision is limited and many ties are present among the flow values.

The OT data or "OTdata" contain flows values over the threshold $u = 2500\,\mathrm{m}$ for the $65$ years period 1913-1977. The historical data or "MAXdata" is simply the $r=12$ largest flows for the period of $143$ years 1770-1912. The exact dates of these events are not known with precision but the years are known and given as comments.

Source

The data were taken from the book by Miquel.

References

Miquel J. (1984) Guide pratique d'estimation des probabilités de crues, Eyrolles (coll. EDF DER).

Parent E. and Bernier J. (2003) Bayesian POT modeling for Historical data. Journal of Hydrology vol. 274, pp. 95-108.

Examples

plot(Garonne)

Description

Translate a (named) vector of GEV parameters into a renewal model.

Usage

gev2Ren(parGev,
        threshold = NULL,
        lambda = NULL,
        w = 1,
        distname.y = c("gpd", "GPD", "lomax", "maxlo"),
        vcovGev = NULL,
        jacobian = TRUE,
        plot = FALSE)

Arguments

Value

A vector of parameters similar to the coefficient vector of an object with class "Renouv". This vector has an element named "lambda" corresponding to the rate of the Homogeneous Poisson Process. The other elements are the parameters for the distribution of POT excesses.

The result has attributes named "distname.y" and "threshold" which can be used to build a Renouv object using the RenouvNoEst function. The result may as well have attributes "jacobian" and "vcov" according to the arguments values. These objects should be used with attention to their element names, since the parameter order may not be the one you expect.

Author(s)

Yves Deville

See Also

The Ren2gev function provide a reciprocal transformation.

Examples

## GEV parameters and block duration
set.seed(1234)
muGev <- rnorm(1); sigmaGev <- rgamma(1, shape = 5)
xiGev <- runif(1, min = -0.2, max = 0.3)
parGev <- c("loc" = muGev, "scale" = sigmaGev, "shape" = xiGev)
parRen <- gev2Ren(parGev, lambda = 1, jacobian = TRUE, plot = TRUE)
## check by reverse transform
parGevCheck <- Ren2gev(parRen, threshold = attr(parRen, "threshold"))
rbind(parGev, parGevCheck)

##=======================================================================
## slightly positive shape convert to "gpd" and "lomax" distributions
##=======================================================================
x <- rgev(n = 100, loc = 3500, scale = 1000, shape = 0.1)
fit.gev <- fgev(x, start = list("loc" = 3000, "scale" = 1000, "shape" = 0.1))
distNames <- c("gpd", "lomax")
namesRen <- c("lambda", "scale", "shape") # works for the 2 target dists
fitNew <- list()
opar <- par(mfrow = c(3, 1))
for (nm in distNames) {  
  parRen <- gev2Ren(parGev = fit.gev$estimate, threshold = 2800,
                    vcov = fit.gev$var.cov, distname.y = nm)
  namesRen <- c("lambda", "scale", "shape")
  myVcov <- attr(parRen, "vcov")[namesRen, namesRen]
  fitNew[[nm]] <- RenouvNoEst(threshold = attr(parRen, "threshold"),
                              estimate = parRen,
                              distname.y = attr(parRen, "distname.y"),
                              cov = myVcov)
  plot(fitNew[[nm]], Tlim = c(1, 200))
}
plot(fit.gev, which = 4)
par(opar)
##=======================================================================
## slightly negative shape convert to "gpd" and "maxlo" distribution
##=======================================================================
x <- rgev(n = 100, loc = 3500, scale = 1000, shape = -0.2)
fit.gev <- fgev(x, start = list("loc" = 3000, "scale" = 1000, "shape" = 0.1))
distNames <- c("gpd", "maxlo")
namesRen <- c("lambda", "scale", "shape") # works for the 2 target dists
fitNew <- list()
opar <- par(mfrow = c(3, 1))
for (nm in distNames) {
  parRen <- gev2Ren(parGev = fit.gev$estimate, threshold = 2800,
                    vcov = fit.gev$var.cov, distname.y = nm)
  myVcov <- attr(parRen, "vcov")[namesRen, namesRen]
  fitNew[[nm]] <- RenouvNoEst(threshold = attr(parRen, "threshold"),
                              estimate = parRen,
                              distname.y = attr(parRen, "distname.y"),
                              cov = myVcov)
  plot(fitNew[[nm]], Tlim = c(1, 200))
}
plot(fit.gev, which = 4)
par(opar)

Description

Goodness-of-fit diagnostics for the distribution of event dates in a (assumed) Poisson process

Usage

gof.date(date,
             start = NULL,
             end = NULL,
             plot = TRUE,
             main = NULL,
             skip = NULL,
             plot.type = "skip")

Arguments

Details

In the homogeneous Poisson process, events occur on a time interval in a uniform fashion. More precisely, for a given time interval the distribution of the event dates conditional to their number $n$ is the distribution of the order statistics of a sample of size $n$ of the uniform distribution on this interval.

When the interval has limits taken at events the uniformity statement remains true, but for inner events. This behaviour is met when start and end are not given and taken as the first and last events in date.

Value

A list

When the number of events corresponding to the indications of args is 0, the function returns NULL with a warning. When the number of events is less than 6 a warning is shown.

Warning

When skipped periods exist the number of events, duration, rate the global KS test must be computed by omitting the skipped periods in the duration and retaining only valid interevents. The indication given in nevt rate and duration should be used only when no skipped period exist (skip = NULL on input) and replaced by effnevt, effrate and effduration otherwise.

Note

In practical contexts missing periods are often met in the datasets. The diagnostic should therefore be applied on every period with no missing data. Even if the event dates seem reasonably uniform, it is a good idea to check that the rates do not differ significantly over intervals.

When some events are missing and no suitable information is given via the skip argument, the global rate, KS.statistic and KS.pvalue are of little interest. Yet the graph might be instructive.

Author(s)

Yves Deville

See Also

interevt function for the determination of interevents ans subsequent diagnostics.

Examples

## Use "Brest" dataset
## simple plot. Kolmogorov-Smirnov is not useful
gof1 <- gof.date(date = Brest$OTdata$date)

## consider missing periods. Much better!
gof2 <- gof.date(date = Brest$OTdata$date,
         skip = Brest$OTmissing,
         start = Brest$OTinfo$start,
         end = Brest$OTinfo$end)

print(gof2$noskip)

## Second type of graph
gof3 <- gof.date(date = Brest$OTdata$date,
         skip = Brest$OTmissing,
         start = Brest$OTinfo$start,
         end = Brest$OTinfo$end,
         plot.type = "omit")

## non-skipped periods at Brest
ns <- skip2noskip(skip = Brest$OTmissing,
                 start = Brest$OTinfo$start,
                 end = Brest$OTinfo$end)

## say 9 plots/diagnostics
oldpar <- par(mar = c(3, 4, 3, 2), mfcol = c(3, 3))

for (i in 1:9) {
  GOF <- gof.date(date = Brest$OTdata$date,
           start = ns$start[i],
           end = ns$end[i])
}

par(oldpar)

Description

Bartlett's goodness-of-fit test for exponential distribution

Usage

gofExp.test(x)

Arguments

Value

A list with elements

Author(s)

Yves Deville

References

See

Yagouti A., Abi-Zeid I., Ouarda, T.B.M.J. and B. Bobée (2001), Revue de processus ponctuels et synthèse de tests statistiques pour le choix d'un type de processus Revue des Sciences de l'Eau, 1, pp. 323-361.

See Also

Among other goodness-of-fit tests ks.test in the stats package. See expplot for a graphical diagnostic.

Examples

## a sample of size 30
 x <- rexp(30)
 res <- gofExp.test(x)

 ## ns samples: p.values should look as uniform on (0, 1)
 ns <- 100
 xmat <- matrix(rexp(30*ns), nrow = ns, ncol = 30)
 p.values <- apply(xmat, 1, function(x) gofExp.test(x)$p.value)
 plot(sort(p.values), type = "p", pch = 16)

Description

Density function, distribution function, quantile function, random generation, hazard and cumulative hazard functions for the Generalised Pareto Distribution.

Usage

dGPD(x, loc = 0.0, scale = 1.0, shape = 0.0, log = FALSE)
   pGPD(q, loc = 0.0, scale = 1.0, shape = 0.0, lower.tail = TRUE)
   qGPD(p, loc = 0.0, scale = 1.0, shape = 0.0, lower.tail = TRUE)
   rGPD(n, loc = 0.0, scale = 1.0, shape = 0.0)
   hGPD(x, loc = 0.0, scale = 1.0, shape = 0.0)
   HGPD(x, loc = 0.0, scale = 1.0, shape = 0.0)

Arguments

Details

Let $\mu$, $\sigma$ and $\xi$ denote loc, scale and shape. The distribution values $y$ are $\mu \leq y < y_{\textrm{max}}$.

When $\xi \neq 0$, the survival function value for $y \geq \mu$ is given by

$S(y) = \left[1 + \xi(y - \mu)/\sigma\right]^{-1/ \xi} \qquad \mu < y < y_{\textrm{max}}$

where the upper end-point is $y_{\textrm{max}} = \infty$ for $\xi >0$ and $y_{\textrm{max}} = \mu -\sigma/ \xi$ for $\xi <0$.

When $\xi = 0$, the distribution is exponential with survival

$S(y) = \exp\left[- (y - \mu)/\sigma\right] \qquad \mu \leq y.$

Value

dGPD gives the density function, pGPD gives the distribution function, qGPD gives the quantile function, and rGPD generates random deviates. The functions hGPD and HGPD return the hazard rate and the cumulative hazard.

Note

The functions are slight adaptations of the [r,d,p,q]gpd functions in the evd package. The main difference is that these functions return NaN when shape is negative, as it might be needed in unconstrained optimisation. The quantile function can be used with p=0 and p=1, then returning the lower and upper end-point.

See Also

fGPD to fit such a distribution by Maximum Likelihood.

Examples

qGPD(p = c(0, 1), shape = -0.2)
shape <- -0.3
xlim <- qGPD(p = c(0, 1), shape = shape)
x <- seq(from = xlim[1], to = xlim[2], length.out = 100)
h <- hGPD(x, shape = shape)
plot(x, h, type = "o", main = "hazard rate for shape < 0")
shape <- 0.2
xlim <- qGPD(p = c(0, 1 - 1e-5), shape = shape)
x <- seq(from = xlim[1], to = xlim[2], length.out = 100)
h <- hGPD(x, shape = shape)
plot(x, h, type = "o", main = "hazard rate shape > 0 ")

Description

Translate a (named) vector of Gumbel parameters into a vector of parameters for a renewal model.

Usage

gumbel2Ren(parGumbel,
           threshold = NULL,
           lambda = NULL,
           w = 1,
           distname.y = c("exponential"),
           vcovGumbel = NULL,
           jacobian = TRUE,
           plot = FALSE)

Arguments

Details

Given a vector of Gumbel parameters and a block duration, there exits an infinity of Renouv models with exponential excesses leading to the prescribed Gumbel distributions for the maximum of the marks on a block with duration w. One of these models may be chosen by specifying either a threshold or a rate lambda.

Value

A vector of Renouv parameters, which can be used with RenouvNoEst.

Caution

All Renouv models lead to the same return level curve whatever be the choice of threshold or lambda. However, when a covariance matrix is given, the covariance matrix for the Renouv parameters and consequently the confidence bounds depend on the threshold or rate. So the computed covariance matrix must in general be considered as putative.

Author(s)

Yves Deville

See Also

gev2Ren for a similar translation from the GEV distribution.

Description

Plotting positions for exponential return level plots.

Usage

Hpoints(n)

Arguments

Details

The plotting positions are numeric values to use as the abscissae corresponding to the order statistics in an exponential return level plot. They range from 1 to about $\log n$. They can be related to the plotting positions given by ppoints.

The returned vector $\mathbf{H}$ has elements

$H_{i} = \frac{1}{n} + \frac{1}{n-1} + \dots + \frac{1}{n + 1 -i}$

for $1 \leq i \leq n$. This is the expectation of the $i$-th order statistic for a sample of the standard exponential distribution, see e.g. chap. 4 of Embrechts et al.

Value

Numeric vector of plotting positions with length n.

Note

For $n$ large enough, the largest value $H_n$ is approximately $\gamma + \log n$ where $\gamma$ is the Euler-Mascheroni constant, and $\exp H_n$ is about $1.78 n$. Thus if the Hpoints are used as plotting positions on a return level plot, the largest observation has a return period of about $1.78 n$ years.

Author(s)

Yves Deville

References

Embrechts P., Klüppelberg C. and Mikosch T. (1997) Modelling Extremal Events for Insurance and Finance. Springer.

See Also

ppoints.

Examples

n <- 30
set.seed(1234)
x <- rGPD(n, shape = 0.2)
plot(exp(Hpoints(n)), sort(x), log = "x",
     main = "Basic return level plot")

Description

Compute a simple (preliminary) estimation for the tree parameters of the mixture of two exponential distributions

Usage

ini.mixexp2(x, plot = FALSE)

Arguments

Details

This function gives estimators using several methods if necessary. The goal is to find the rates rate1, rate2 and the mixing probability prob1 with the 'feasibility' constraints 0 < rate1 < rate2 and 0 < prob1 < 1.

First the method of moments is used. If the estimates are feasible they are returned with method = "moments". If not, the estimates are derived using two linear regressions. A regression without constant using only the smallest values gives an estimator of the mean rate. A regression using only the largest values gives rate1 and prob1. Yet the constraints must be fulfilled. If they are, the estimates are returned (together with method = "Hreg" suggesting a cumulative hazard regression). If not, a (poor) default estimate is returned with method = "arbitrary".

Value

A list

Note

The method of moments is implemented in mom.mixexp2. Further investigations are needed to compare the estimators (moments or Hreg) and select the best strategy.

Note that this function returns the estimate within a list and no longer as a vector with named elements as was the case before.

Author(s)

Yves Deville

See Also

See MixExp2, mom.mixexp2.

Examples

set.seed(1234)
x <- rmixexp2(n = 100, prob1 = 0.5, rate2 = 4)
res <- ini.mixexp2(x, plot = TRUE)

Description

Compute intervent durations from events dates

Usage

interevt(date,
            skip = NULL, noskip = NULL)

Arguments

Details

Interevents are the time differences between successive dates. When the date argument contains occurrence times $T_i$ for successive events of an homogeneous Poisson process, interevents $T_i -T_{i-1}$ are mutually independent with the same exponential distribution.

When some time intervals are skipped independently from the event point process, we may consider the interevents $T_i-T_{i-1}$ between two non-skipped events such that the time interval $(T_{i-1},\,T_i)$ does not contains any skipped interval. These interevents still are mutually independent with the same exponential distribution. When skip or noskip is not NULL the computation therefore only retains couples of two successive datetimes "falling" in the same non-skipped period, which number can therefore be associated with the interevent.

Value

A list mainly containing a interevt data.frame.

Note

Only one of the two arguments skip and noskip should be given in the call. In each case, the rows of the returned data.frame objects describe periods in chronological order. That is: start at row 2 must be after the end value of row 1 and so on.

Note that there are usually less interevents than dates since two successive dates will be retained for an interevent only when they are not separated by missing period. As a limit case, there can be no interevents if the noskip periods contain only one date from the date vector.

Author(s)

Yves Deville

See Also

gof.date for goodness-of-fit diagnostics for dates of events expplot for diagnostics concerning the exponential distribution.

Examples

## Use Brest data
ie <- interevt(date = Brest$OTdata$date, skip = Brest$OTmissing)

expplot(ie$interevt$duration, rate = 1 / mean(ie$interevt$duration),
  main = "No threshold")

## keep only data over a threshold
ind1 <- Brest$OTdata$Surge >= 35
ie1 <- interevt(Brest$OTdata$date[ind1], skip = Brest$OTmissing)
expplot(ie1$interevt$duration, main = "Threshold = 35")

## increase threshold
ind2 <- Brest$OTdata$Surge >= 55
ie2 <- interevt(date = Brest$OTdata$date[ind2], skip = Brest$OTmissing)
expplot(ie2$interevt$duration, main = "Threshold = 55 cm")

Description

Jackson's statistic for the exponentiality test.

Usage

Jackson(x, norm = FALSE)

Arguments

Details

The value(s) of the statistic are the ratio of two weighted means of the order statistics.

Value

A numeric vector of length 1 when x is a vector, or with length nrow(x) when x is a matrix.

References

J. Beirlant and T. de Weit and Y. Goegebeur(2006) A Goodness-of-fit Statistic for Pareto-Type Behaviour, J. Comp. Appl. Math., 186(1), pp. 99-116

Description

Jackson's test of exponentiality

Usage

Jackson.test(x, method = c("num", "sim", "asymp"), nSamp = 15000)

Arguments

Details

Compute the Jackson's test of exponentiality. The test statistic is the ratio of weighted sums of the order statistics. Both sums can also be written as weighted sums of the scalings.

The Jackson's statistic for a sample of size $n$ of the exponential distribution can be shown to be approximately normal. More precisely $\sqrt{n}(J_n -2)$ has approximately a standard normal distribution. This distribution is used to compute the $p$-value when method is "asymp". When method is "num", a numerical approximation of the distribution is used. Finally, when method is "sim" the $p$-value is computed by simulating nSamp samples of size length(x) and estimating the probability to have a Jackson's statistic larger than that of the 'observed' x.

Value

A list of results.

Note

Jackson's test of exponentiality works fine for a Lomax alternative (GPD with heavy tail). It then reaches nearly the same power as a Likelihood Ratio (LR) test, see Kozubowski et al. It can be implemented more easily than the LR test because simulated values of the test statistic can be obtained quickly enough to compute the $p$-value by simulation.

Author(s)

Yves Deville

References

J. Beirlant and T. de Weit and Y. Goegebeur(2006) "A Goodness-of-fit Statistic for Pareto-Type Behaviour", J. Comp. Appl. Math., 186(1), pp. 99-116.

T.J. Kozubowski, A. K. Panorska, F. Qeadan, A. Gershunov and D. Rominger (2009) "Testing Exponentiality Versus Pareto Distribution via Likelihood Ratio" Comm. Statist. Simulation Comput. 38(1), pp. 118-139.

See Also

The Jackson function computing the statistic and the LRExp.test function.

Examples

set.seed(1234)
x <- rGPD(n = 50, loc = 0, scale = 1, shape = 0.1)
Jackson.test(x, method = "num")$p.value
Jackson.test(x, method = "asymp")$p.value
Jackson.test(x, method = "sim")$p.value

Description

Log-likelihood, AIC, BIC and number of observations of an object of class "Renouv".

Usage

## S3 method for class 'Renouv'
AIC(object, ..., k = 2)
## S3 method for class 'Renouv'
BIC(object, ...)
## S3 method for class 'Renouv'
logLik(object, ...)
## S3 method for class 'Renouv'
nobs(object, ...)

Arguments

Caution

Comparing log-likelihoods, AIC or BIC for different Renouv objects makes sense only when these share the same data and the same threshold.

Note

logLik, AIC and BIC can be used with an object of class "Renouv" which makes use of historical data. In this case, the number of observations may be misleading since a single historical observation may concern dozens of years and thus have a much greater impact on the estimation of the tail than an "ordinary" observation.

Author(s)

Yves Deville

See Also

The AIC, nobs generic functions.

Description

Density function, distribution function, quantile function and random generation for the Lomax distribution.

Usage

dlomax(x, scale = 1.0, shape = 4.0, log = FALSE)
   plomax(q, scale = 1.0, shape = 4.0, lower.tail = TRUE)
   qlomax(p, scale = 1.0, shape = 4.0)
   rlomax(n, scale = 1.0, shape = 4.0)

Arguments

Details

The Lomax distribution function with shape $\alpha > 0$ and scale $\beta > 0$ has survival function

$S(y) = \left[1 + y/\beta \right]^{-\alpha} \qquad (y > 0)$

This distribution has increasing hazard and decreasing mean residual life (MRL). The coefficient of variation decreases with $\alpha$, and tends to $1$ for large $\alpha$. The default value $\alpha=4$ corresponds to $\textrm{CV} = \sqrt{2}$.

Value

dlomax gives the density function, plomax gives the distribution function, qlomax gives the quantile function, and rlomax generates random deviates.

Note

This distribution is sometimes called log-exponential. It is a special case of Generalised Pareto Distribution (GPD) with positive shape $\xi > 0$, scale $\sigma$ and location $\mu=0$. The Lomax and GPD parameters are related according to

$\alpha = 1/\xi, \qquad \beta = \sigma/\xi.$

The Lomax distribution can be used in POT to describe excesses following GPD with shape $\xi>0$ thus with decreasing hazard and increasing Mean Residual Life.

Note that the exponential distribution with rate $\nu$ is the limit of a Lomax distribution having large scale $\beta$ and large shape $\alpha$, with the constraint on the shape/scale ratio $\alpha/\beta = \nu$.

References

Johnson N. Kotz S. and N. Balakrishnan Continuous Univariate Distributions vol. 1, Wiley 1994.

Lomax distribution in Wikipedia

See Also

flomax to fit the Lomax distribution by Maximum Likelihood.

Examples

shape <- 5; scale <- 10
xl <- qlomax(c(0.00, 0.99), scale = scale, shape = shape)
x <- seq(from = xl[1], to = xl[2], length.out = 200)
f <- dlomax(x, scale = scale, shape = shape)
plot(x, f, type = "l", main = "Lomax density")
F <- plomax(x, scale = scale, shape = shape)
plot(x, F, type ="l", main ="Lomax distribution function")

Description

Likelihood Ratio statistic for the exponential distribution vs. GPD.

Usage

LRExp(x, alternative = c("lomax", "GPD", "gpd", "maxlo"))

Arguments

Details

The Likelihood-Ratio statistic is actually $W:=-2 \log \textrm{LR}$ where LR is the ratio of the likelihoods exponential to alternative distribution.

Value

The LR statistic value.

Note

When the alternative is "lomax" or "maxlo", the statistic has a distribution of mixed type under the null hypothesis of exponentiality. This is a mixture of a distribution of continuous type (with positive values) and of a Dirac mass at LR = 0. The probability mass $\textrm{Pr}\{\textrm{LR} = 0\}$ can be computed using the pGreenwood1 function. More precisely, the probability mass is pGreenwood1(n) for the Lomax alternative and 1.0 - pGreenwood1(n) for the maxlo alternative, where n is the sample size length(x).

See Also

LRExp.test for the related LR test of exponentiality.

Description

Likelihood Ratio test of exponentiality vs. GPD.

Usage

LRExp.test(x,
              alternative = c("lomax", "GPD", "gpd", "maxlo"),
              method = c("num", "sim", "asymp"),
              nSamp = 15000,
              simW = FALSE)

Arguments

Details

The Lomax and maxlo alternatives correspond to a GPD alternative with positive shape parameter $\xi > 0$ (Lomax) and GPD with $\xi < 0$ (maxlo).

The asymptotic distribution of the Likelihood-ratio statistic is known. For the GPD alternative, this is a chi-square distribution with one df. For the Lomax alternative, this is the distribution of a product $BC$ where $B$ and $C$ are two independent random variables following a Bernoulli distribution with probability parameter $p = 0.5$ and a chi-square distribution with one df.

• When method is "num", a numerical approximation of the distribution is used. This method is not unlike that used by Kozubowski et al., but a different approximation is used. However, if x has a length $n > 500$, the method is turned to "asymp".

• When method is "sim", nSamp samples of the exponential distribution with the same size as x are drawn and the LR statistic is computed for each sample. The $p$-value is simply the estimated probability that a simulated LR is greater than the observed LR.

• Finally when method is "asymp", the asymptotic distribution is used.

Value

A list of results with elements statistic, p.value and method. Other elements are

Note

For the Lomax alternative, the distribution of the test statistic has mixed type: it can take any positive value as well as the value $0$ with a positive probability mass. The probability mass is the probability that the ML estimate of the GPD shape parameter is negative, and a good approximation of it is provided by the pGreenwood1 function. Note that this probability converges to its limit $0.5$ very slowly, which suggests that the asymptotic distribution provides poor results for medium sample sizes, say $< 100$.

Similarly for a maxlo alternative, the distribution of the test statistic has mixed type: it can take any positive value as well as the value $0$ with a positive probability mass approximately given by 1 -pGreenwood1(n) where $n$ is the sample size.

Author(s)

Yves Deville

References

T.J. Kozubowski, A. K. Panorska, F. Qeadan, A. Gershunov and D. Rominger (2009) "Testing Exponentiality Versus Pareto Distribution via Likelihood Ratio" Comm. Statist. Simulation Comput. 38(1), pp. 118-139.

The approximation method used is described in the Renext Computing Details report.

See Also

Lomax, Maxlo, GPD for the alternatives used here.

Examples

set.seed(1234)
x <- rGPD(n = 50, loc = 0, scale = 1, shape = 0.1)
LRExp.test(x, method = "num")$p.value
LRExp.test(x, method = "asymp")$p.value
## Not run: 
## requires much time
LRExp.test(x, method = "sim")$p.value

## End(Not run)

Description

Likelihood Ratio statistic for the Gumbel distribution vs. GEV.

Usage

LRGumbel(x, alternative = c("frechet", "GEV"))

Arguments

Details

The Likelihood-Ratio statistic is actually $W:=-2 \log \textrm{LR}$ where LR is the ratio of the likelihoods Gumbel to alternative distribution.

Value

The LR statistic value.

Note

When the alternative is "frechet", the statistic has a distribution of mixed type under the null hypothesis of a Gumbel distribution.

Author(s)

Yves Deville

See Also

LRGumbel.test for the related LR test of Gumbelity.

Description

Likelihood Ratio test of Gumbel vs. GEV

Usage

LRGumbel.test(x,
                 alternative = c("frechet", "GEV"),
                 method = c("num", "sim", "asymp"),
                 nSamp = 1500,
                 simW = FALSE)

Arguments

Details

The asymptotic distribution of the Likelihood-ratio statistic is known. For the GEV alternative, this is a chi-square distribution with one df. For the Fréchet alternative, this is the distribution of a product $XY$ where $X$ and $Y$ are two independent random variables following a Bernoulli distribution with probability parameter $p = 0.5$ and a chi-square distribution with one df.

• When method is "num", a numerical approximation of the distribution is used.

• When method is "sim", nSamp samples of the Gumbel distribution with the same size as x are drawn and the LR statistic is computed for each sample. The $p$-value is simply the estimated probability that a simulated LR is greater than the observed LR. This method requires more computation time than the tow others.

• Finally when method is "asymp", the asymptotic distribution is used.

Value

A list of results with elements statistic, p.value and method. Other elements are

Note

For the Fréchet alternative, the distribution of the test statistic has mixed type: it can take any positive value as well as the value $0$ with a positive probability mass. The probability mass is the probability that the ML estimate of the GEV shape parameter is negative.

When method is "sim", the computation can be slow because each of the nSamp simulated values requires two optimisations. The "asymp" method provides an acceptable precision for $n \geq 50$, and may even be used for $n \geq 30$.

Author(s)

Yves Deville

Examples

set.seed(1234)
x <- rgumbel(60)
res <- LRGumbel.test(x)

Description

Density function, distribution function, quantile function and random generation for the 'maxlo' distribution.

Usage

dmaxlo(x, scale = 1.0, shape = 4.0, log = FALSE)
   pmaxlo(q, scale = 1.0, shape = 4.0, lower.tail = TRUE)
   qmaxlo(p, scale = 1.0, shape = 4.0)
   rmaxlo(n, scale = 1.0, shape = 4.0)

Arguments

Details

The 'maxlo' distribution function with shape $\alpha>0$ and scale $\beta>0$ is a special case of Generalised Pareto (GPD) with negative shape $\xi < 0$ and location at zero. This is the finite upper endpoint case of the GPD. Its name is nonstandard and was chosen to suggest some form of symmetry with respect to the Lomax distribution.

The survival function is

$S(y) = \left[1-y/\beta\right]^\alpha \qquad 0 < y < \beta$

This distribution has a coefficient of variation smaller than $1$.

Value

dmaxlo gives the density function, pmaxlo gives the distribution function, qmaxlo gives the quantile function, and rmaxlo generates random deviates.

Note

The 'maxlo' and GPD parameters are related according to

$\alpha = -1/\xi, \qquad \beta = -\sigma/\xi.$

where $\sigma$ is the scale parameter of the GPD. Since only GPD with $\xi > -0.5$ seem to be used in practice, this distribution should be used with $\alpha > 2$.

This distribution can be used in POT to describe bounded excesses following GPD with shape $\xi < 0$. The scale parameter $\beta$ then represents the upper end-point of the excesses, implying the finite upper end-point $u + \beta$ for the levels, where $u$ is the threshold. It can be used in Renouv with a fixed scale parameter, thus allowing a control of the upper end-point.

This distribution is simply a rescaled version of a beta distribution and also a rescaled version of a Kumaraswamy distribution. The name "maxlo" is used here to suggest a form of symmetry to Lomax distribution.

See Also

fmaxlo to fit such a distribution by Maximum Likelihood.

Examples

xs <- rmaxlo(500, shape = 2.2, scale = 1000)
hist(xs, main = "'maxlo' distribution"); rug(xs)

xs <- rmaxlo(500, shape = 4, scale = 1000)
hist(xs, main = "'maxlo' distribution"); rug(xs)

x <- seq(from = -10, to = 1010, by = 2)
plot(x = x, y = dmaxlo(x, shape = 4, scale = 1000),
     type = "l", ylab = "dens",
     col = "orangered", main = "dmaxlo and dgpd")
abline(h = 0)
lines(x = x, y = dgpd(x, shape = -1/4, scale = 250),
     type = "l",
     col = "SpringGreen3", lty = "dashed")

Description

Probability functions associated to the mixture of two exponential distributions.

Usage

dmixexp2(x, prob1,
            rate1 = 1.0, rate2 = rate1 + delta, delta,
            log = FALSE)
   pmixexp2(q, prob1,
            rate1 = 1.0, rate2 = rate1 + delta, delta,
            log = FALSE)
   qmixexp2(p, prob1,
            rate1 = 1.0, rate2 = rate1 + delta, delta)
   rmixexp2(n, prob1,
            rate1 = 1.0, rate2 = rate1 + delta, delta)
   hmixexp2(x, prob1,
            rate1 = 1.0, rate2 = rate1 + delta, delta)
   Hmixexp2(x, prob1,
            rate1 = 1.0, rate2 = rate1 + delta, delta)

Arguments

Details

The density function is the mixture of two exponential densities

$f(x) = \alpha_1 \lambda_1 \, e^{-\lambda_1 x} + (1-\alpha_1) \lambda_2 \, e^{-\lambda_2x} \qquad x > 0$

where $\alpha_1$ is the probability given in prob1 while $\lambda_1$ and$\lambda_2$ are the two rates given in rate1 and rate2.

A 'naive' identifiability constraint is

$\lambda_1 < \lambda_2$

i.e. rate1 < rate2, corresponding to the simple constraint delta > 0. The parameter delta can be given instead of rate2.

The mixture distribution has a decreasing hazard, increasing Mean Residual Life (MRL) and has a thicker tail than the usual exponential. However the hazard, MRL have a finite non zero limit and the distribution behaves as an exponential for large return levels/periods.

The quantile function is not available in closed form and is computed using a dedicated numerical method.

Value

dmiwexp2, pmiwexp2, qmiwexp2, evaluates the density, the distribution and the quantile functions. dmixexp2 generates a vector of n random draws from the distribution. hmixep2 gives hazard rate and Hmixexp2 gives cumulative hazard.

Examples

rate1 <- 1.0
rate2 <- 4.0
prob1 <- 0.8
qs <- qmixexp2(p = c(0.99, 0.999), prob1 = prob1,
               rate1 = rate1, rate2 = rate2) 
x <- seq(from = 0, to = qs[2], length.out = 200)
F <- pmixexp2(x, prob1 = prob1, rate1 = rate1, rate2 = rate2)
plot(x, F, type = "l", col = "orangered", lwd = 2,
     main = "Mixexp2 distribution and quantile for p = 0.99")
abline(v = qs[1])
abline(h = 0.99)

Description

Compute the moment estimation for the tree parameters of the mixture of two exponential distributions

Usage

mom.mixexp2(x)

Arguments

Details

The three parameters (probability and the two rates) are computed from the first three moments (theoretical and sample). It can be shown that the inverse rates are obtained solving a quadratic equation. However the roots can be negative or complex and the estimates are not valid ones.

Value

A list with elements