cPTdensity: Nonparametric Bayesian Estimation of Joint Scedasis Density...

Description Usage Arguments Details Value Author(s) References Examples

View source: R/cPTdensity.R

Description

This function computes a Bayesian joint scedasis density estimate generating a posterior sample of a finite Mixture of Polya trees.

Usage

1
2
cPTdensity(XY, tau = 0.95, raw = TRUE, prior, mcmc, state,status, data =
          sys.frame(sys.parent()), na.action = na.fail)

Arguments

XY

data frame from which the estimate is to be computed; first column corresponds to time, the second and third columns correspond to the variables of interest X and Y respectively.

tau

value used to threshold the data z = min(X, Y); by default threshold = quantile(z, tau).

raw

logical; if TRUE, X and Y will be converted to unit Fréchet scale. If FALSE, X and Y will be understood as already in unit Fréchet scale. If a single variable is provided raw is automatically set to TRUE.

prior

a list giving the prior information. The list includes the following parameter: a0 and b0 giving the hyperparameters for prior distribution of the precision parameter of the Poly tree prior, alpha giving the value of the precision parameter (it must be specified if alpha is missing, see details below), optionally M giving the finite level to be considered (if M is specified, a partially specified mixture of Polya trees model is fitted), tau1 and tau2 giving the hyperparameters of the inverted gamma prior distribution for the centering Beta scale parameter,m0 and S0 giving the hyperparameters of the inverted gamma prior distribution for the centering Beta scale parameter, and al, be giving the value of the shape and scale parameter of the centering distribution.

mcmc

a list giving the MCMC parameters. The list must include the following integers: nburn giving the number of burn-in scans, nskip giving the thinning interval, nsave giving the total number of scans to be saved, ndisplay giving the number of saved scans to be displayed on screen (the function reports on the screen when every ndisplay iterations have been carried out), tune1, tune2, and tune3 giving the positive Metropolis tuning parameter for the baseline shape, scale, and precision parameter, respectively (the default value is 1.1)

state

a list giving the current value of the parameters. This list is used if the current analysis is the continuation of a previous analysis.

status

a logical variable indicating whether this run is new (TRUE) or the continuation of a previous analysis (FALSE). In the latter case the current value of the parameters must be specified in the object state.

data

data frame.

na.action

a function that indicates what should happen when the data contain NAs. The default action (na.fail) causes PTdensity to print an error message and terminate if there are any incomplete observations.

Details

This function learns about the joint scedasis function using a Mixture of Polya Trees (MPT) prior as proposed in Palacios and de Carvalho (2020). In the particular case where XY contains no third column, the function learns about the scedasis function of Einmahl et al (2016) using an MPT prior. The details are as follows. Let

Z_i = min(X_i, Y_i)

. The model is given by:

Z1,...,Zn | G ~ G

G | alpha,a,b ~ PT(Pi^{al,be},\textit{A})

where, the the PT is centered around a Beta(al,be) distribution, by taking each m level of the partition Π^{al, be} to coincide with the k/2^m, k=0,…,2^m quantile of the Beta(al,be) distribution. The family \textit{A}={alphae: e \in E^{*}}, where E^{*}=\bigcup_{m=0}^{M} E^m and E^m is the m-fold product of E=\{0,1\}, was specified as alpha{e1 … em}=α m^2.

In the univariate case,the following proper priors can be assigned:

al | m0, S0 ~ LNormal(m0,S0)

be | tau1, tau2 ~ LNormal(tau1,tau2)

To complete the model specification, independent hyperpriors are assumed,

alpha | a0, b0 ~ Γ(a0,b0)

The precision parameter, alpha, of the PT prior can be considered as random, having a gamma distribution, Γ(a0,b0), or fixed at some particular value. To let alpha to be fixed at a particular value, set a0 to NULL in the prior specification.

In the computational implementation of the model, Metropolis–Hastings steps are used to sample the posterior distribution of the baseline and precision parameters.

Value

c

Joint scedasis density estimator.

k

number of exceedances above the threshold.

w

standardized indices of exceedances.

Y

raw data.

al

giving the value of the baseline shape parameter.

be

giving the value of the baseline scale parameter.

alpha

giving the value of the precision parameter.

The plot method depicts the estimated joint scedasis density.

Author(s)

Miguel de Carvalho and Vianey Palacios

References

Palacios, V., de Carvalho, M. (2020) Bayesian semiparametric modeling of jointly heteroscedastic extremes. Preprint.

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
## Not run: 
## Example from Palacios and de Carvalho (2020, submitted)
library(evd)
## Initial state
state <- NULL
## MCMC parameters
nburn <- 2000
nsave <- 1000
nskip <- 0
ndisplay <- 500
mcmc <- list(nburn=nburn,nsave=nsave,nskip=nskip,ndisplay=ndisplay,
 tune1=1.1,tune2=1.1,tune3=1.1)
## Prior information
prior<-list(a0=1,b0=1,M=8,m0=.01,S0=.01,tau1=.01,tau2=.01);

T <- 5000
time <- seq(1/T, 1, by = 1/T)
set.seed(12333)
aux <- matrix(0, T, 2)
for (i in 1:T) {
    aux[i,]<-rbvevd(1, dep =.5, asy=c(sin(pi*time[i]),time[i]),
    model="alog",mar1 = c(1, 1, 1), mar2 = c(1, 1, 1))
}
XY <- cbind(time, aux)
fit <- cPTdensity(XY, prior = prior, mcmc = mcmc, state = state, status =
TRUE)
plot(fit)

## End(Not run)

extremis documentation built on Nov. 27, 2020, 9:07 a.m.