pFailure | R Documentation |
This function computes the empirical estimate of the probability of falling into two types of failure regions.
pFailure(n, beta, u1, u2, mar1, mar2, type, plot, xlab, ylab,
nlevels=10)
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 ( |
plot |
A logical value; if |
xlab , ylab |
A character string equivalent representing the graphical parameters as in |
nlevels |
The number of contour levels to be displayed. |
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\}
A list containing AND
and/or OR
, depending on the type
argument. Each element is a length(u1) x length(u2)
matrix.
Simone Padoan, simone.padoan@unibocconi.it, https://faculty.unibocconi.it/simonepadoan/; Boris Beranger, borisberanger@gmail.com https://www.borisberanger.com;
Marcon, G., Naveau, P. and Padoan, S.A. (2017) A semi-parametric stochastic generator for bivariate extreme events Stat, 6, 184-201.
dExtDep
, rExtDep
, fExtDep
, fExtDep.np
# 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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.