gibbs_statio: Gibbs Sampler for GEV (MCMC)

In proto4426/PissoortThesis: Thesis in Extreme Value Theory

Usage

gibbs_mcmc.own(start, nbr.chain = length(start), propsd = c(0.4, 0.1, 0.1),
  Gumbel = F, iter = 2000, burnin = ceiling(iter/2 + 1),
  data = max_years$data)

gibbs_mcmc.own_WithoutGumbel(start, nbr.chain = length(start),
  propsd = c(0.4, 0.1, 0.1), iter = 2000, burnin = ceiling(iter/2 + 1),
  data = max_years$data)

Arguments

start	numeric vector of length 3 containing the starting values for the parameters theta= (location, LOG-scale and shape). It is advised explore different ones, and typically take the MPLE
iter	The number of iterations of the algorithm. Must e high enough to ensure convergence
burnin	Determines value for burn-in
data	numeric vector containing the GEV in block-maxima
proposd	The proposal's standard deviations : controlling the cceptance rate. To facilitate convergence, it is advised to target an acceptance rate of around 0.25 when all components of theta are updated simultaneously, and 0.40 when the components are updated one at a time. is as proposed by Coles (2001) but we should tune this value. (from package ismev)

Value

A named list containing

n.chains : The number of chains generated melted in a data.frame

mean.acc_rates : the meanS of the acceptance rates

out.chain : The generated chainS

dic.vals : contains the DIC values (for further diagnostics on predictive accuracy, see ?dic)

out.ind : The generated individual chainS (in a list)

Author(s)

Antoine Pissoort, [e-mail]

Examples

data("max_years")
fn <- function(par, data) -log_post0(par[1], par[2], par[3], data)
param <- c(mean(max_years$df$Max),log(sd(max_years$df$Max)), 0.1 )
# opt <- optim(param, fn, data = max_years$data,
#              method="BFGS", hessian = TRUE)
opt <- nlm(fn, param, data = max_years$data,
           hessian=T, iterlim = 1e5)
start <- opt$estimate
Sig <- solve(opt$hessian)
ev <- eigen( (2.4/sqrt(2))^2 * Sig)
varmat <- ev$vectors %*% diag(sqrt(ev$values)) %*% t(ev$vectors)
## GIBBS
 set.seed(100)
iter <- 2000
gibb1 <- gibbs_mcmc.own(start, iter = iter)

proto4426/PissoortThesis documentation built on May 26, 2019, 10:31 a.m.