Learn R Programming

ExtremalDep (version 0.0.4-2)

pFailure: 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)

Value

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

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\} $$

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

Run this code

# 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")

# \donttest{
  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)
# }

Run the code above in your browser using DataLab