gibbs.Nstatio1: Gibbs Sampler for nonstationary GEV (MCMC)
In proto4426/PissoortThesis: Thesis in Extreme Value Theory

Description Usage Arguments Value Author(s) Examples

Compute the Gibbs sampler accounting for nonstationarity (trend) in GEV. It is computed based on a diffuse normal prior.

log_post1(mu0, mu1, logsig, xi, data, model.mu = mu0 + mu1 * tt,
  mnpr = c(30, 0, 0, 0), sdpr = c(40, 40, 10, 10), rescale.time = T,
  ic = F)

gibbs.trend.own(start, propsd = c(0.5, 2.5, 0.08, 0.08), iter = 1000,
  burnin = ceiling(iter/2 + 1), data = max_years$data,
  keep.same.seed = NULL, rescale.time = T, Progress.Shiny = NULL,
  .mnpr = c(30, 0, 0, 0), .sdpr = c(40, 40, 10, 10))

`data`	numeric vector containing the GEV in block-maxima
`start`	named list of length 4 (this number determines the number of chains generated) containing the starting values for the parameters theta=(intercept mu0, trend mu1, 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
`keep.same.seed`	sets a seed at each iterations that is the integer you specify times the iteration number.
`Progress.Shiny`	Additional parameter that handles the progress bar in the pertaining Shiny application.
`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. It must be wel chosen (e.g. Trial-and-error method)

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)

Antoine Pissoort, antoine.pissoort@student.uclouvain.be

data("max_years")
data <- max_years$data

fn <- function(par, data) -log_post1(par[1], par[2], par[3],
                                     par[4], data)
param <- c(mean(max_years$df$Max), 0, log(sd(max_years$df$Max)), -0.1 )
opt <- optim(param, fn, data = max_years$data,
             method = "BFGS", hessian = T)

# Starting Values
set.seed(100)
start <- list() ; k <- 1
while(k < 5) { # starting value is randomly selected from a distribution
  # that is overdispersed relative to the target
  sv <- as.numeric(rmvnorm(1, opt$par, 50 * solve(opt$hessian)))
  svlp <- log_post1(sv[1], sv[2], sv[3], sv[4], max_years$data)
  print(svlp)
  if(is.finite(svlp)) {
    start[[k]] <- sv
    k <- k + 1
  }
}

# k chains with k different starting values
set.seed(100)
gibbs.trend <- gibbs.trend.own(start, propsd = c(.5, 1.9, .15, .12),
                               iter = 1000)
colMeans(do.call(rbind, gibbs.trend$mean_acc.rates))

param.chain <- gibbs.trend$out.chain[ ,1:4]

### Plot of the chains
chains.plotOwn(gibbs.trend$out.chain )