returns {ExtremalDep} | R Documentation |
Predicts the probability of future simultaneous exceedances
returns(mcmc, summary.mcmc, y, plot=FALSE)
mcmc |
The output of the |
summary.mcmc |
The output of the |
y |
A 2-column matrix of unobserved thresholds. |
plot |
If |
Computes for a range of unobserved extremes (larger than those observed in a sample), the pointwise mean from the posterior predictive distribution of such predictive values. The probabilities are calculated through
P(Y_1 > y_1, Y_2 > y_2) = \frac{2}{k} \sum_{j=0}^{k-2} (\eta_{j+1} - \eta_j)
\times \left(
\frac{(j+1) B(y_1/(y_1+y_2)| j+2, k-j-1)}{y_1}
- \frac{(k-j-1) B(y_2/(y_1+y_2)| k-j, j+1)}{y_2}
\right),
where B(x|a,b)
denotes the cumulative distribution function of a Beta random variable with shape a,b>0
. See Marcon et al. (2016, p.3323) for details.
Returns a vector whose length is equal to the number of rows of the input value y
.
Simone Padoan, simone.padoan@unibocconi.it, https://mypage.unibocconi.it/simonepadoan/; Boris Beranger, borisberanger@gmail.com https://www.borisberanger.com/; Giulia Marcon, giuliamarcongm@gmail.com
Marcon G., Padoan, S.A. and Antoniano-Villalobos I. (2016) Bayesian Inference for the Extremal Dependence. Electronic Journal of Statistics, 10.2, 3310-3337.
if (interactive()){
# This reproduces some of the results shown in Fig. 1 (Marcon, 2016).
set.seed(1890)
data <- evd::rbvevd(n=100, dep=.6, asy=c(0.8,0.3), model="alog", mar1=c(1,1,1))
nsim = 500000
burn = 400000
mu.nbinom = 3.2
var.nbinom = 4.48
hyperparam <- list(a.unif=0, b.unif=.5, mu.nbinom=mu.nbinom, var.nbinom=var.nbinom)
k0 = 5
pm0 = list(p0=0.06573614, p1=0.3752118)
eta0 = ExtremalDep:::rcoef(k0, pm0)
mcmc <- bbeed(data, pm0, eta0, k0, hyperparam, nsim,
prior.k = "nbinom", prior.pm = "unif")
w <- seq(0.001, .999, length=100)
summary.mcmc <- summary.bbeed(w, mcmc, burn, nsim, plot=TRUE)
plot.bbeed(type = "A", x=w, mcmc=mcmc, summary.mcmc, nsim=nsim, burn=burn)
plot.bbeed(type = "h", x=w, mcmc=mcmc, summary.mcmc, nsim=nsim, burn=burn)
plot.bbeed(type = "pm", x=w, mcmc=mcmc, summary.mcmc, nsim=nsim, burn=burn)
plot.bbeed(type = "k", x=w, mcmc=mcmc, summary.mcmc, nsim=nsim, burn=burn)
y <- seq(10,100,2)
y <- as.matrix(expand.grid(y,y))
probs <- returns(mcmc = mcmc, summary.mcmc = summary.mcmc, y = y, plot = TRUE)
}