Package {mevr}

Version:	1.3.3
Date:	2026-10-13
Title:	Fitting the Metastatistical Extreme Value Distribution MEVD
Maintainer:	Harald Schellander <harald.schellander@geosphere.at>
Depends:	R (≥ 4.1), dplyr, rlang
Encoding:	UTF-8
Imports:	graphics, stats, EnvStats, parallel, foreach, doParallel, bamlss, mgcv, mmand, data.table (≥ 1.18.0),
Description:	Extreme value analysis with the metastatistical extreme value distribution MEVD (Marani and Ignaccolo, 2015, <doi:10.1016/j.advwatres.2015.03.001>) and some of its variants. In particular, analysis can be performed with the simplified metastatistical extreme value distribution SMEV (Marra et al., 2019, <doi:10.1016/j.advwatres.2019.04.002>) and the temporal metastatistical extreme value distribution TMEV (Falkensteiner et al., 2023, <doi:10.1016/j.wace.2023.100601>). Parameters can be estimated with probability weighted moments, maximum likelihood and least squares. The data can also be left-censored prior to a fit. Density, distribution function, quantile function and random generation for the MEVD, SMEV and TMEV are included. In addition, functions for the calculation of return levels including confidence intervals are provided. For a description of use cases please see the provided references.
License:	GPL-3
Suggests:	testthat (≥ 3.0.0)
Config/testthat/edition:	3
URL:	https://haraldschellander.github.io/mevr/
Config/roxygen2/version:	8.1.1
NeedsCompilation:	no
Packaged:	2026-10-13 09:58:23 UTC; hschellander
Author:	Harald Schellander [aut, cre] (Package creator and maintainer), Alexander Lieb [ctb] (Coded first versions of MEVD and SMEV), Marc-Andre Falkensteiner [ctb] (Developed the TMEV)
Repository:	CRAN
Date/Publication:	2026-10-13 10:20:02 UTC

Fit Weibull distribution to censored data

Description

Finds the optimal left-censoring threshold(s) at which the data series should be censored to make sure that the observations in the tail are likely sampled from a Weibull distribution

Usage

censored_weibull_fit(x, thresholds, warn)

Arguments

Value

A tibble with the optimal threshold itself and the Weibull scale and shape parameters obtained from the censored sample.

Examples

data("dailyrainfall")
wbtest <- weibull_tail_test(dailyrainfall)
censored_weibull_fit(wbtest, 0.9)

Daily rainfall data

Description

A dataset containing daily rainfall intended to be used with the package mevr

Usage

data(dailyrainfall)

Format

The dataset contains real world daily rainfall observations from a station in the northern Alps. The series contains values from 1971 to 1985 and are assumed to be Weibull distributed. This data series is intended to be used as is as input data for the package mevr to fit the metastatistical extreme value distribution and its variants with different estimation methods.

The dataset is a dataframe with two columns, dates and val:

dates

Days of class Date in the format YYYY-MM-DD

val

Rainfall observations corresponding to the date in the row. The value is the 24 hour sum from the morning hours of day-1 to the morning hours of day.

Examples

## Load example data
data(dailyrainfall)

## explore dataset
head(dailyrainfall)
hist(dailyrainfall$val)
plot(dailyrainfall$val, type = "o")

The Metastatistical Extreme Value Distribution

Description

Density, distribution function, quantile function and random generation for the MEV distribution with shape parameter 'w', scale parameter 'c'. Parameter 'n' refers either to the mean number of wet days per year in case of the simplified MEV, or to the number of wet days for each year.

Usage

dmev(x, w, c, n)

pmev(q, w, c, n)

qmev(p, w, c, n)

rmev(N, w, c, n)

Arguments

Value

dmev gives the density function, pmev gives the distribution function, qmev gives the quantile function and rmev provides random realizations of the SMEV and MEVD.

Functions

• pmev(): distribution quantile function

• qmev(): quantile function

• rmev(): random generation function

Examples

# SMEV
dmev(1200:1300, 0.7, 20, 110)
pmev(1200:1300, 0.7, 70, 110)
qmev(1 - 1 / seq(5,50,5), 0.7, 70, 110)

# MEVD: 50-year event of 10 years Weibull series
w <- rnorm(10, 0.8, 0.1) # shape parameter
c <- rnorm(10, 200, 30) # scale parameter
n <- rnorm(10, 200, 50) # number of wet days
qmev(1 - 1 / 50, w, c, n)

# rl-plot
rp <- seq(5, 50, 5)
rl <- qmev(1 - 1 / rp, w, c, n)
pp <- (1:length(rp)) / (length(rp) + 1)
plot(pp, rl, type = "o")

The non-stationary Metastatistical Extreme Value Distribution

Description

Quantile function for the TMEV distribution with a Weibull parent distribution.

Usage

dtmev(x, data)

ptmev(q, data)

qtmev(p, data)

Arguments

Value

dtmev gives the density function, ptmev gives the distribution function, and qtmev gives the quantile function of the TMEV.

Functions

• ptmev(): distribution quantile function

• qtmev(): distribution quantile function

Examples

data(dailyrainfall)
fit <- ftmev(dailyrainfall)
rp <- pp.weibull(fit$maxima)
rl <- qtmev(1 - 1 / rp, fit$data)
plot(rp,
     sort(fit$maxima),
     main = "TMEV",
     ylab = "return level",
     xlab = "return period (years)")
lines(rp, rl, type = "l")

Detect Rainfall Events with Specified Separation Time (Dry Time) and Minimum Rain Threshold

Description

This function identifies rainfall events in time-series data based on a minimum rainfall threshold and a separation time between events. It uses morphological operations to detect and separate events, and it can ignore short events that do not meet a minimum duration criterion. Based on Marra et al. 2024: "Predicting Extreme Sub-Hourly Precipitation Intensification Based on Temperature Shifts“. Hydrology and Earth System Sciences 28, Nr. 2: 375–89. doi:10.5194/hess-28-375-2024.

Usage

event_separation(
  data,
  separation_period = 360,
  time_resolution = 10,
  ignore_event_duration = 30,
  min_rain = 0.1
)

Arguments

Details

The function works by first dilating the rainfall data using a structuring element defined by the separation time, followed by erosion to remove noise. It then detects events based on the separation time and the minimum rainfall threshold, returning the time intervals where events occurred. Events that span more than two years or are shorter than the 'ignore_event_duration' are ignored.

Value

A list with the following elements:

• 'data': A modified version of the input data with an additional column 'is_event', indicating whether each time step is part of a detected event (1 for event, 0 for no event). This is useful in the next step ordinary_events.

• 'fromto': A tibble with two columns, 'from' and 'to', representing the indices of the start ('from') and end ('to') of each detected rainfall event.

• 'ts_res': The time resolution used in the function (in minutes).

Examples

## Not run: 
# Example data
mock_data <- data.frame(
  groupvar = seq.POSIXt(Sys.time(),
                        by = "10 min",
                        length.out = 100),
  val = c(rep(0, 10), runif(20, 0, 1), rep(0, 70))
)

# Detect rainfall events with a separation of 6 hours
# and a minimum rain threshold of 0.1
result <- event_separation(mock_data,
                           separation_period = 360,
                           time_resolution = 10,
                           min_rain = 0.1)
print(result$fromto)  # View detected events

## End(Not run)

Fitting the Metastatistical Extreme Value Distribution (MEVD)

Description

Fit the MEVD distribution to rainfall observations with different estimation methods.

Usage

fmev(
  data,
  threshold = 0,
  method = c("pwm", "mle", "ls"),
  censor = FALSE,
  censor_opts = list(),
  warn = TRUE,
  left_cens = 0
)

Arguments

Details

With the aim of weakening the requirement of an asymptotic assumption for the GEV distribution, a metastatistical approach was proposed by Marani and Ignaccolo (2015). The MEVD is defined in terms of the distribution of the statistical parameters describing "ordinary" daily rainfall occurrence and intensity. The MEVD accounts for the random process of event occurrence in each block and the possibly changing probability distribution of event magnitudes across different blocks, by recognizing the number of events in each block, n, and the values of the shape and scale parameters w and C of the parent Weibull distribution to be realisations of stochastic variables. The MEVD can then be written as

\begin{aligned} F = \frac{1}{T} \sum_{j=1}^T \prod_{k \in A_j} &\left( 1 - \exp\left( -\left(\frac{x}{C(j,k)}\right)^{w(j,k)} \right) \right) \end{aligned}

for w > 0 and C > 0. With T fully recorded years, yearly C and w can be estimated by fitting a Weibull distribution to the values x of this year, and n is the number of ordinary events per year, i.e. all rainfall events larger than a threshold.

If the probability distribution of daily rainfall is assumed to be time-invariant, the MEVD can be simplified to

F = [1 - exp(-x/C)^w]^n

with single values for the shape and scale parameters w and C. n is then the mean number of wet days at this location (Marra et al., 2019; Schellander et al., 2019).

As is shown e.g. Schellander et al., 2019, probability weighted moments should be preferred over maximum likelihood for the estimation of the Weibull parameters w and C. Therefore method = 'pwm' is the default.

The MEVD can also be used for sub-daily precipitation (Marra et al., 2019). In that case n has to be adapted accordingly to the 'mean number of wet events' per year.

This function returns the parameters of the fitted MEVD distribution as well as some additional fitting results and input parameters useful for further analysis.

Value

A list of class mevr with components:

c

Vector of Weibull scale parameters of the MEVD, where each component refers to one year.

w

Vector of Weibull shape parameters of the MEVD, where each component refers to one year.

n

Number of wet events per year. Wet events are defined as rainfall > threshold.

params

A named vector of the fitted parameters.

maxima

Maximum values corresponding to each year.

data

Filtered data used to fit the MEVD along with additional components useful for further analysis.

years

Vector of years as YYYY.

threshold

The chosen threshold.

method

Method used to fit the MEVD.

type

The type of distribution ("MEVD").

Author(s)

Harald Schellander, Alexander Lieb

References

Marani, M. and Ignaccolo, M. (2015) 'A metastatistical approach to rainfall extremes', Advances in Water Resources. Elsevier Ltd, 79(Supplement C), pp. 121-126. doi:10.1016/j.advwatres.2015.03.001.

Schellander, H., Lieb, A. and Hell, T. (2019) 'Error Structure of Metastatistical and Generalized Extreme Value Distributions for Modeling Extreme Rainfall in Austria', Earth and Space Science, 6, pp. 1616-1632. doi:10.1029/2019ea000557.

See Also

fsmev, ftmev

Examples

data(dailyrainfall)
fit <- fmev(dailyrainfall, method = "mle")
fit
plot(fit)

Fitting the simplified Metastatistical Extreme Value Distribution (SMEV)

Description

Fit the SMEV distribution to rainfall observations with different estimation methods.

Usage

fsmev(
  data,
  threshold = 0,
  method = c("pwm", "mle", "ls"),
  censor = FALSE,
  censor_opts = list(),
  warn = TRUE,
  sd = FALSE,
  sd.method = "boot",
  R = 502,
  left_cens = 0
)

Arguments

Details

The SMEV was introduced by (Marra et al., 2019) as a simplified version of the MEVD (Marani and Ignaccolo, 2015) with the assumption of a stationary parent Weibull distribution as

F = [1 - exp(-x/C)^w]^n

for w > 0 and C > 0 being the Weibull shape and scale parameter and n > 0 being the mean number of wet days over all years. Wet days are defined as rainfall events > threshold. As it was shown by e.g. Schellander et al., 2019, probability weighted moments should be preferred over maximum likelihood for the estimation of the Weibull parameters w and C. Therefore method = 'pwm' is the default.

Confidence intervals of the SMEV distribution can be calculated using a non parametric bootstrap technique. Note that this very slow.

The estimators above assume that all ordinary events are Weibull distributed. This is usually acceptable for daily rainfall, but not for sub-daily durations, where the lower part of the sample is affected by the gauge resolution and by weak events that do not share the tail of the intense ones. Fitted to the complete sample, the shape parameter is then overestimated and return levels are strongly underestimated. Following Marra et al. (2019, 2020) the Weibull distribution can be fitted to the tail only by left-censoring: with left_cens = 0.9 the smallest 90 percent of the ordinary events are censored, i.e. their magnitudes are ignored but they keep their non-exceedance probability, and n is still the mean number of all ordinary events per year. With method = 'ls' this is the least squares regression in Weibull-transformed coordinates used by Marra et al. (2020), method = 'mle' maximizes the likelihood of the left-censored sample and method = 'pwm' uses partial probability weighted moments (Wang, 1990). For sub-daily data, ordinary events should be defined with event_separation and ordinary_events.

This function returns the parameters of the fitted SMEV distribution as well as some additional fitting results and input parameters useful for further analysis.

Value

A list of class mevr with components:

c

Single value of the Weibull scale parameter of the SMEV.

w

Single value of the Weibull shape parameter of the SMEV.

n

Mean number of wet events, averaged over all years. Wet events are defined as rainfall > threshold.

params

A named vector of the fitted parameters.

maxima

Maximum values corresponding to each year.

std

Standard error of fitted parameters (if sd = TRUE).

varcov

Covariance matrix of fitted parameters (if sd = TRUE).

data

Filtered data used to fit the SMEV and additional components useful for further analysis.

years

Vector of years as YYYY.

threshold

The chosen threshold.

method

Method used to fit the MEVD.

censor

TRUE when the data series was left-censored and FALSE otherwise.

left_cens

The fixed left-censoring quantile used for estimating c and w.

type

The type of distribution ("SMEV").

rejected

If censor = TRUE, TRUE when the Weibull tail assumption is rejected, and FALSE otherwise. Omitted if censor = FALSE.

Author(s)

Harald Schellander, Alexander Lieb

References

Marra, F. et al. (2019) 'A simplified MEV formulation to model extremes emerging from multiple nonstationary underlying processes', Advances in Water Resources. Elsevier Ltd, 127(April), pp. 280-290. doi:10.1016/j.advwatres.2019.04.002.

Marra, F. et al. (2020) 'A unified framework for extreme subdaily precipitation frequency analyses based on ordinary events', Geophysical Research Letters, 47, e2020GL090209. doi:10.1029/2020GL090209.

Wang, Q. J. (1990) 'Estimation of the GEV distribution from censored samples by method of partial probability weighted moments', Journal of Hydrology, 120, pp. 103-114.

See Also

fmev, ftmev

Examples

data(dailyrainfall)

fit <- fsmev(dailyrainfall)
fit
plot(fit)

# estimate c and w from the largest 25 percent of the ordinary events only
fit_lc <- fsmev(dailyrainfall, left_cens = 0.75)
fit_lc

# left censor data prior to fitting
set.seed(123)
sample_dates <- seq.Date(from = as.Date("2000-01-01"),
                         to = as.Date("2020-12-31"),
                         by = 1)
sample_data <- data.frame(dates = sample_dates,
                          val = sample(rnorm(length(sample_dates))))
d <- sample_data |>
  filter(val >= 0 & !is.na(val))

fit <- fsmev(d)
fit_c <- fsmev(d,
               censor = TRUE,
               censor_opts = list(thresholds = c(seq(0.5, 0.9, 0.1), 0.95),
                                  mon = 1,
                                  nrtrials = 2,
                                  R = 100))

rp <- 2:100
rl <- return.levels.mev(fit, return.periods = rp)
rl_c <- return.levels.mev(fit_c, return.periods = rp)
plot(sort(pp.weibull(fit$maxima)), sort(fit$maxima),
  xlab = "Return period (a)", ylab = "daily rain (mm)")
lines(rl$rp, rl$rl)
lines(rl_c$rp, rl_c$rl, col = "red")
legend("bottomright", legend = c("std", "censored"),
  col = c("black", "red"), lty = c(1, 2), lwd = c(1, 1.5), bty = "n")

Fitting the temporal Metastatistical Extreme Value Distribution (TMEV)

Description

Fit the temporal MEVD distribution TMEV to rainfall observations with a cyclic spline to account for seasonality.

Usage

ftmev(
  data,
  threshold = 0,
  minyears = 10,
  day_year_interaction = FALSE,
  verbose = FALSE,
  yday_ti_shape_k = 10,
  yday_ti_scale_k = 10,
  year_ti_shape_k = 10,
  year_ti_scale_k = 10
)

Arguments

Details

With the aim of exploiting the full temporal information for parameter estimation, Falkensteiner et al., (2023) introduced the TMEV, which is an explicitly non-stationary formulation of the MEVD (Marani and Ignaccolo, 2015). Adopting a Weibull distribution for ordinary rainfall events, the assumption of yearly constant coefficients is relaxed by allowing the Weibull parameters to fluctuate with time. The TMEV can then be written as

\begin{aligned} F = \frac{1}{T} \sum_{j=1}^T &\prod_{k \in A_j} \left( 1 - \exp\left( -\left(\frac{x}{C(j,k)}\right)^{w(j,k)} \right) \right) \end{aligned}

with w > 0 and C > 0 being the Weibull shape and scale parameters, and $A_j \subseteq (1, ..., 366)$ being the wet days in year j. The temporal and the superimposed seasonal dependence on w and c is modeled with a cyclic seasonal effect on the day of the year.

Technically this is accomplished by fitting a cyclic spline to the daily rainfall values. The following formula is used for the fitting procedure of both the Weibull scale and shape parameter with the function bamlss from package bamlss:

parameter = x \sim s(year) + ti(yday, bs = "cc", k = 10)

The first effect models the long-term temporal trend of the parameter with a thin-plate spline. The second effect models the superimposed seasonal fluctuations of the parameter with the 'day of the year' with a cyclic cubic regression spline and 10 knots, to ensure a smooth transition between December and January. The number of knots (k) in the above equation can be set separately for the year and yday effect as well as separately for the shape and scale parameter of the Weibull distribution. This can be done by overwriting the parameters yday_ti_shape_k, yday_ti_scale_k, year_ti_shape_k, year_ti_scale_k in the call to ftmev. Note that these values depend on many factors, such as the structure of the data, the TMEV is fitted to. Please refer to the documentation of the packages bamlss and, in particular, mgcv.

For data series with lengths < 10 years, the first temporal effect is changed to a simple linear time trend.

For trend analysis, an additional interaction term can be added to the model formula. The following term models the relationship between the seasonality as day of the year and the year itself with a combination of a thin plate and a cyclic cubic spline:

\begin{aligned} ti(&year, yday, bs = c("tp", "cc"), \\ &d = c(1, 1), k = c(year\_ti\_k, yday\_ti\_k)) \end{aligned}

This function returns the parameters of the fitted TMEV distribution as well as some additional fitting results and input parameters useful for further analysis.

Value

A list of class mevr with components:

c

Vector of Weibull scale parameters of the TMEV, where each row refers to an event (daily rainfall).

w

Vector of Weibull shape parameters of the TMEV, where each row refers to an event (daily rainfall).

n

Number of wet events per year. Wet events are defined as rainfall > threshold.

maxima

Maximum values corresponding to each year.

data

A data frame with the data used to fit the TMEV and the fitted Weibull parameters c and w.

years

Vector of years as YYYY.

threshold

The chosen threshold.

x

The fitted bamlss object.

type

The type of distribution ("TMEV").

minyears

The minimum number of years used to fit the TMEV as provided.

Author(s)

Marc-Andre Falkensteiner, Harald Schellander

References

Marani, M. and Ignaccolo, M. (2015) 'A metastatistical approach to rainfall extremes', Advances in Water Resources. Elsevier Ltd, 79(Supplement C), pp. 121-126. doi:10.1016/j.advwatres.2015.03.001.

Falkensteiner, M., Schellander, H., Hell, T. (2023) 'Accounting for seasonality in the metastatistical extreme value distribution', (Weather and Climate Extremes, 42, 2023, doi:10.1016/j.wace.2023.100601).

See Also

fmev, fsmev

Examples

data(dailyrainfall)
filtered_data <- subset(
  dailyrainfall,
  as.numeric(format(dates, "%Y")) >= 1971 &
  as.numeric(format(dates, "%Y")) <= 1978
)
fit <- ftmev(filtered_data, minyears = 5)
plot(fit, type = "rl")

# temporal trend of the Weibull parameters
pred <- predict(fit)
pred_year <- predict(fit, term = "year")
boxplot(c.pred ~ year, data = pred)
with(pred_year, lines(year - 1970,
                      c.pred.year,
                      type = "b",
                      pch = 20,
                      col = "red"))

Identifies ordinary rainfall events by calculating the maximum within rainfall events defined with event_separation

Description

This function calculates the rolling sum of values in the provided data over a specified duration and identifies the events with the highest sum within the defined time intervals. It is useful for detecting significant events based on aggregated values over time.

Usage

ordinary_events(x, duration, na.rm = TRUE)

Arguments

Value

Returns a tibble with individual rainfall events that can be used as input for functions fsmev, fmev, ftmev.

Examples

## Not run: 
# Example usage
x <- list(
  data = data.frame(val = runif(100),
  groupvar = seq.POSIXt(Sys.time(),
                        by = "10 min",
                        length.out = 100)
                        ),
  ts_res = 10,
  fromto = data.frame(from = c(1, 51), to = c(50, 100))
)
duration <- 30
result <- ordinary_events(x, duration)
print(result)

## End(Not run)

Plot graphs of MEVD, SMEV or TMEV fit

Description

A return level plot, qq-plot, pp-plot and a histogram with the fitted density is produced

Usage

## S3 method for class 'mevr'
plot(
  x,
  q = c(2, 10, 20, 30, 50, 75, 100, 150, 200),
  ci = FALSE,
  type = c("all", "rl", "qq", "pp", "hist"),
  ...
)

Arguments

Value

No return value, only a plot is produced.

Examples

data(dailyrainfall)

# fit a simplified MEVD
fit <- fsmev(dailyrainfall)
fit
plot(fit)

# fit MEVD
fit <- fmev(dailyrainfall, method = "ls")
fit
plot(fit)

Weibull plotting position

Description

Calculates the Weibull plotting position for the given maxima

Usage

pp.weibull(x)

Arguments

Value

A numeric vector of Weibull plotting positions corresponding to the given maxima x

Examples

data(dailyrainfall)
fit <- fsmev(dailyrainfall)
rp <- pp.weibull(fit$maxima)
rl <- return.levels.mev(fit, return.periods = rp)
plot(rp,
     sort(fit$maxima),
     xlab = "Return period (years)",
     ylab = "Return level",
     main = fit$type)
lines(rp, rl$rl)

TMEV prediction

Description

Takes a mevr object where the TMEV has been fitted to rainfall data and calculates bamlss predictions for the distributional parameters and the model terms. Basically a wrapper to the corresponding function predict.bamlss

Usage

## S3 method for class 'mevr'
predict(object, newdata, term, ...)

Arguments

Details

See also the details of ftmev for an explanation of the model terms used to fit the temporal trend of the Weibull parameters. The basis dimensions yday_ti_shape_k, yday_ti_scale_k, year_ti_shape_k, year_ti_scale_k are taken from the fitting process, i.e. the call to ftmev.

Value

A data.frame with the supplied covariables and the predicted parameters.

See Also

ftmev

Examples

data(dailyrainfall)

# restrict for the sake of speed
idx <- which(as.POSIXlt(dailyrainfall$date)$year + 1900 < 1976)
data <- dailyrainfall[idx, ]

f <- ftmev(data, minyears = 5)
predict(f, term = "year")

Print method for object of class mevr

Description

Print nicely formatted output of the fit to the MEVD and its variants

Usage

## S3 method for class 'mevr'
print(x, digits = max(3, getOption("digits") - 3), ...)

Arguments

Value

A nicely formatted output of the fitting results.

Examples

data(dailyrainfall)

# fit a simplified MEVD
fit <- fsmev(dailyrainfall)
print(fit)

Return Levels for the MEVD/SMEV/TMEV extreme value distributions

Description

Calculate return levels for a MEVD, SMEV or TMEV extreme value distributions from an object of class mevr.

Usage

return.levels.mev(
  x,
  return.periods = c(2, 10, 20, 30, 50, 75, 100, 150, 200),
  ci = FALSE,
  alpha = 0.05,
  method = "boot",
  R = 502,
  ncores = 2L
)

Arguments

Details

Note that bootstraping the confidence intervals is very slow.

Value

A list with return levels, chosen return periods and, if ci=TRUE, alpha/2 and 1 - alpha/2 confidence intervals.

Examples

data(dailyrainfall)

fit <- fmev(dailyrainfall)
return.levels.mev(fit)
plot(fit)

Weibull tail test

Description

This functions provides a way to test if observed rainfall maxima from a data series are likely samples from a parent distribution with a Weibull tail. The concept and the code is based on the paper 'Marra F, W Amponsah, SM Papalexiou, 2023. Non-asymptotic Weibull tails explain the statistics of extreme daily precipitation. Adv. Water Resour., 173, 104388, doi:10.1016/j.advwatres.2023.104388. They also provide the corresponding Matlab code (https://zenodo.org/records/7234708).

Usage

weibull_tail_test(
  data,
  threshold = 0,
  mon = 1,
  cens_quant = 0.9,
  p_test = 0.1,
  R = 500
)

Arguments

Details

Null-Hyothesis: block maxima are samples from a parent distribution with Weibull tail (tail defined by a given left-censoring threshold). If the fraction of observed block maxima outside of the interval defined by p_test is larger than p_test the null hypothesis is rejected.

Value

A tibble with the test outcome and other useful results:

Examples

data("dailyrainfall")
weibull_tail_test(dailyrainfall)

# generate data
set.seed(123)
sample_dates <- seq.Date(from = as.Date("2000-01-01"),
                         to = as.Date("2010-12-31"),
                         by = 1)
sample_data <- data.frame(dates = sample_dates,
                          val = sample(rnorm(length(sample_dates))))
d <- sample_data |>
  filter(val >= 0 & !is.na(val))
fit_uncensored <- fsmev(d)

# censor the data
thresholds <- c(seq(0.1, 0.9, 0.1), 0.95)
p_test <- 0.1
res <- lapply(thresholds, function(x) {
  weibull_tail_test(d, cens_quant = x, p_test = p_test, R = 200)
})
res <- do.call(rbind, res)

# find the optimal left-censoring threshold
cens <- censored_weibull_fit(res, thresholds)
cens$optimal_threshold
cens$quantile

# plot return levels censored vs uncensored
rp <- c(2:100)
rl_uncensored <- return.levels.mev(fit_uncensored, return.periods = rp)$rl
rl_censored <- qmev(1 - 1/rp, cens$shape, cens$scale, fit_uncensored$n)
plot(rp,
     rl_uncensored,
     type = "l",
     log = "x",
     ylim = c(0, max(rl_censored, rl_uncensored)),
     ylab = "return level",
     xlab = "return period (a)")
points(pp.weibull(fit_uncensored$maxima), sort(fit_uncensored$maxima))
lines(rp, rl_censored, type = "l", col = "red")
legend("bottomright", legend = c("uncensored",
        paste0("censored at ", round(cens$optimal_threshold, 1), "mm")),
        col = c("black", "red"), lty = c(1, 1))