extrememix (version 0.0.1)

fmgpd: MGPD Estimation

Description

Fit of the MGPD model using an MCMC algorithm.

Usage

fmgpd(x, it, k, start = NULL, var = NULL, prior = NULL, thin = 1, burn = 0)

Value

fmgpd returns a list with three elements:

• chain: a matrix of size (it - burn)/thin\(\times\)5, reporting in each column the posterior sample for each parameter.

• var: a matrix of size it\(\times\)5 reporting the variance of the proposal distribution for each parameter.

• data: the dataset used for estimation.

Arguments

x

A vector of positive observations.

it

Number of iterations of the algorithm.

k

number of mixture components for the bulk. Must be either 2, 3, or 4.

start

A list of starting parameter values.

var

A list of starting proposal variance.

prior

A list of hyperparameters for the prior distribution.

thin

Thinning interval.

burn

Burn-in.

Details

Estimation of the MGPD is carried out using an adaptive block Metropolis-Hastings algorithm. As standard, the user needs to specify the data to use during estimation, the number of mixture components for the bulk, the number of iterations of the algorithm, the burn-in period (by default equal to zero) and the thinning interval (by default equal to one). To run the algorithm it is also needed the choice of the starting values, the starting values of the proposal variances, and the parameters of the prior distribution. If not provided, these are automatically set as follows:

• starting values: \(u\) is chosen by the function ithresh of the threshr package; \(\xi\) and \(\sigma\) are chosen by the fpot function of evd for data over the threshold; \(\mu\) and \(\eta\) are chosen as estimates of the gammamixEM function from the mixtools package; \(w\) is chosen as the vector with entries \(1/k\).

• variances: variances for \(\sigma\) and \(u\) are chosen as the standard deviation of the normal distribution whose 0.01 quantile is equal to 0.9 times the starting value of the associated parameter. The parameters \(\mu_i\) and \(\eta_i\) are sampled jointly and the proposal variance is chosen using the same method as for \(\sigma\) with respect to the parameter \(\mu\). The proposal variance for \(w\) is 0.1 and the proposal variance for \(\xi\) is 0.1 if negative and 0.25 if positive.

• prior distributions: the prior distribution for \(\xi\) and \(\sigma\) is the objective prior $$p(\xi,\sigma) = \sigma^{-1}(1+\xi)^{-1}(1+2\xi)^{-1/2}.$$ The prior for the threshold \(u\) is Normal with mean chosen as for the starting values and the standard deviation is chosen so that the 0.05 quantile of the prior is equal to the median of the data. The priors for the parameters \(\mu_i\) and \(\eta_i\) are Gammas with mean chosen as for the starting values and shapes equal to 0.001 so to give high prior variances. The prior for the weigths is the non-informative Dirichlet with parameter 1.

The user can also select any of the three inputs above.

• The starting values start must be a list with entries xi, sigma, u, mu, eta and w. The length of mu, eta and w must be k.

• The variances var must be a list with entries xi, sigma, u, mu and w. The length of mu must be k.

• The prior prior must be a list with entries u, mu_mu, mu_eta, eta_mu and eta_eta. u gives the mean and the standard deviation of the Normal prior for \(u\). The vectors of length k, mu_mu and eta_mu give the prior means of \(\mu\) and \(\eta\), whilst mu_eta and eta_eta give their prior shapes.

References

Behrens, Cibele N., Hedibert F. Lopes, and Dani Gamerman. "Bayesian analysis of extreme events with threshold estimation." Statistical Modelling 4.3 (2004): 227-244.

do Nascimento, Fernando Ferraz, Dani Gamerman, and Hedibert Freitas Lopes. "A semiparametric Bayesian approach to extreme value estimation." Statistics and Computing 22.2 (2012): 661-675.

See Also

fggpd, mgpd

Examples

# \donttest{
data(rainfall)
## Small number of iterations and burn-in for quick execution
model1 <- fmgpd(rainfall, k = 2, it = 250, burn = 50, thin = 25)
start <- list(xi = 0.2, sigma = 2, u = 10, mu = c(2,5), eta = c(2,2) , w = c(0.4,0.6))
var <- list(xi = 0.01, sigma = 1, u = 3, mu = c(3,3), w = 0.01)
prior <- list(u = c(22,5), mu_mu = c(2,5), mu_eta = c(0.01,0.01),
         eta_mu = c(3,3),eta_eta = c(0.01,0.01))

model2 <- fmgpd(rainfall, k= 2, it = 250, start = start, var =var, prior = prior)
# }