pFailure {ExtremalDep}

Probability of falling into a failure region

Description

This function computes the empirical estimate of the probability of falling into two types of failure regions.

Usage

pFailure(n, beta, u1, u2, mar1, mar2, type, plot, xlab, ylab, 
         nlevels=10)

Arguments

n

An integer indicating the number of points generated to compute the empirical probability.

beta

A vector representing the Bernstein polynomial coefficients.

u1, u2

Vectors of positive reals representing some high thresholds.

mar1, mar2

Vectors of marginal (GEV) parameters

type

A character string indicating if the failure region includes at least one exceedance ("or"), or both exceednaces ("and"). If "both" then probabilities to fall in both regions are computed.

plot

A logical value; if TRUE contour plots of the probalities are displayed.

xlab, ylab

A character string equivalent representing the graphical parameters as in par.

nlevels

The number of contour levels to be displayed.

Details

The two failure regions are:

A_u = \left\{ (v_1, v_2): v_1>u_1 \textrm{ or } v_2 > u_2\right\}

and

B_u = \left\{ (v_1, v_2): v_1>u_1 \textrm{ and } v_2 > u_2\right\}

where (v_1,v_2) \in R_+^2 and u_1, u_2>0.

Exceedances samples y_{1,1}, \ldots, y_{n_1} and y_{1,2}, \ldots, y_{n_2} are generating according to Algorithm 3 of Marcon et al. (2017) and the the estimates of the probability of falling in A_u and B_u are given by

\hat{P}_{A_u} = \frac{1}{n}\sum{i=1}^n I\left\{ y_{i,1}> u_1^* \textrm{ or } y_{i,2}> u_2^* \right\}

and

\hat{P}_{B_u} = \frac{1}{n}\sum{i=1}^n I\left\{ y_{i,1}> u_1^* \textrm{ and } y_{i,2}> u_2^* \right\}

Value

A list containing AND and/or OR, depending on the type argument. Each element is a length(u1) x length(u2) matrix.

Author(s)

Simone Padoan, [e-mail], https://faculty.unibocconi.it/simonepadoan/; Boris Beranger, [e-mail] https://www.borisberanger.com;

References

Marcon, G., Naveau, P. and Padoan, S.A. (2017) A semi-parametric stochastic generator for bivariate extreme events Stat, 6, 184-201.

See Also

dExtDep, rExtDep, fExtDep, fExtDep.np

Examples

# Example wind speed and wind gust

data(WindSpeedGust)
years <- format(ParcayMeslay$time, format="%Y")
attach(ParcayMeslay[which(years %in% c(2004:2013)),])
  
WS_th <- quantile(WS,.9)
DP_th <- quantile(DP,.9)
  
pars.WS <- evd::fpot(WS, WS_th, model="pp")$estimate
pars.DP <- evd::fpot(DP, DP_th, model="pp")$estimate
  
data_uf <- trans2UFrechet(cbind(WS,DP), type="Empirical")
  
rdata <- rowSums(data_uf)
r0 <- quantile(rdata, probs=.90)
extdata <- data_uf[rdata>=r0,]
  
SP_mle <- fExtDep.np(method="Frequentist", data=extdata, k0=10, 
                     type="maxima")


  pF <- pFailure(n=50000, beta=SP_mle$Ahat$beta,
                 u1=seq(from=19, to=28, length=200), mar1=pars.WS,
                 u2=seq(from=40, to=60, length=200), mar2=pars.DP, 
                 type="both", plot=TRUE, 
                 xlab="Daily-maximum Wind Speed (m/s)",
                 ylab="Differential of Pressure (mbar)", nlevels=15)