Scripts-R/Bayes_evdbayes.R

In proto4426/PissoortThesis: Thesis in Extreme Value Theory

setwd('/home/piss/Documents/Extreme/R resources/IRM')
load("data1.Rdata")

# Bayesian Analysis constrcuted with functions in package evdbayes
##################################################################
library(evdbayes)
# This package uses the Metropolis algorithm for MCMC

library(PissoortThesis)
library(coda)



mat <- diag(c(10000, 10000, 100))
# Large variance prior : near flat
pn <- prior.norm(mean = c(0,0,0), cov = mat)

n <- 1000 ;   t0 <- c(31, 1 ,0) ;   s <- c(.4,.1,.1)
# s contains sd for proposal distributions
max.mc <- posterior(n, t0, prior = pn, lh = "gev",
                    data = max_years$data, psd = s)
max.mc

mcmc.max <- mcmc (max.mc, start = 0, end = 1000)
plot(mcmc.max, den = F, sm = F)

maxpst <- mposterior(init = t0, prior = pn, lh = "gev", method = "BFGS",
                     data = max_years$data)
round(maxpst$par, 2)
# use this a initial value, then redo the analysis (...)
max.mc <- posterior(n, maxpst$par, prior = pn, lh = "gev",
                    data = max_years$data, psd = s)

maxpsann <- mposterior(init = t0, prior = pn, lh = "gev", method = "SANN",
                       data = max_years$data)

# psd ?

psd <- ar.choice(init = t0, prior = pn, lh = "gev", data = max_years$data,
                 psd = rep(0.1,3), tol = rep(0.02, 3))$psd
max.mc2 <- posterior(n, maxpst$par, prior = pn, lh = "gev",
                     data = max_years$data, psd = psd)
mcmc.maxF <- mcmc (max.mc2, start = 0, end = 1000)
plot(mcmc.maxF, den = F, sm = F)


### Diagnostics (1)

mcmc.diag <- window(mcmc.maxF, start = 200)
gewek_par <- geweke.diag(mcmc.diag)
2*dnorm(gewek_par$z)  # High pvals -> do not RHO that means of both chains are =
geweke.diag(mcmc.diag, .2, .4)    # (explore values)
# Ok, no abs(>2) so burn-in of 200 seems acceptable.
geweke.plot(mcmc.maxF)

raftery.diag(mcmc.diag, r = 0.01, s = 0.85)  # tune the values
#to obtain required accuracy levels. It seems that xi needs much longer chain

autocorr.plot(mcmc.maxF)   ; autocorr(mcmc.maxF)
# try with 5000 values

max.mc3 <- posterior(5000, maxpst$par, prior = pn, lh = "gev",
                     data = max_years$data, psd = psd)
mcmc.maxF2 <- mcmc (max.mc3, start = 0, end = 1000)

bwf <- function(x) sd(x)/2  # From symmetry of the densities ?
plot(mcmc.maxF2, trace = FALSE, bwf = bwf)
summary(mcmc.maxF2)


######## Point Process (with threshold obv) approach

#igamma(shape = c(38.9,7.1,47), scale = c(1.5,6.3,2.6))
priorpp <- prior.quant(shape = c(0.05, 0.25, 0.1), scale = c(1.5,3,2)) # change

n <- 10000  ;  s <- c(2, .35, .07) ;   b <- 2000  # change

summary(rn.prior <- posterior(n, t0, priorpp, "none", psd = s, burn = b))
#t0 <- c(43.2, 7.64, 0.32) ; s <- c(2, .2, .07)
rn.post <- posterior(n, maxpst$par, priorpp, data = max_all, "pp", thresh = 32,
                     noy = 116, psd = s, burn = b)
plot(mc.post <- mcmc(rn.post, start = 0, end = 8000))


### add TREND in location  ( + see later)

pn.t <- prior.norm(mean = c(0,0,0), cov = mat, trendsd = 500)
rn.post.t <- posterior(50000, c(maxpst$par,1), pn.t, data = max_all, "pp", thresh = 32,
                     noy = 116, psd = c(s,.25), burn = b)
summary(rn.post.t)
plot( mc.post.t <- mcmc(rn.post.t, start = 0, end = 8000) )
plot(rn.post.t)

# library(mvtnorm)
# library(KernSmooth)
# tmp <- bkde2D(rn.post.t[,1:2], bandwidth = c(0.34,0.34))
# contour(seq(min(rn.post.t[,1]), max(rn.post.t[,1]), 0.04), seq(min(rn.post.t[,2]), max(rn.post.t[,2]),0.04), tmp$fhat,
#         drawlabels = FALSE, #levels = seq(0,0.5,0.04),
#         #xlim = c(208.18,210.22), ylim = c(17.5,25.0),
#         cex.lab = 1.25, cex.axis = 1.25, lwd = 2,
#         xlab = "Alpha Parameter", ylab = "Beta Parameter")


## Return Leels

rl.pred(max.mc2, npy = 365.25, qlim = c(30,45))
rl.pred(max.mc3, qlim = c(30,45))




######  Reduce dependence in the chain (see large (p)acf) : THNINNING
# Store only k-th value generated by the chain.

thin.post <- posterior(10000, maxpst$par, priorpp, data = max_all, "pp", thresh = 32,
                       noy = 116, psd = c(s, .25), burn = b, thin = 5)
summary(thin.post)
plot(mc.post <- mcmc(thin.post, start = , end = nrow(thin.post)))

thin.post <- posterior(10000, c(maxpst$par,1), pn.t, data = max_all, "pp", thresh = 32,
                       noy = 116, psd = c(s, .25), burn = b, thin = 5)
summary(thin.post)
plot(mc.post <- mcmc(thin.post, start = , end = nrow(thin.post)))



####### Multiple Chains : Diagnostic (2) of convergence : Gelman
# Run chains with different starting values

init1 <- c(30, 2, -.3) ;   init2 <- c(33, 1, 0)  ;  init3 <- c(31, 1.5, -.15)
mc1 <- posterior(5000, init1, prior = pn, lh = "gev",
                 data = max_years$data, psd = psd)
mc2 <- posterior(5000, init2, prior = pn, lh = "gev",
                 data = max_years$data, psd = psd)
mc3 <- posterior(5000, init3, prior = pn, lh = "gev",
                 data = max_years$data, psd = psd)
post.mcl <- mcmc.list(mcmc(mc1), mcmc(mc2), mcmc(mc3))
traceplot(post.mcl)  # mixing is good for all
gelman.diag(post.mcl, transform = TRUE)
gelman.plot(post.mcl)
# Everything seems okay here, for this MCMC.



#### For the QUANTILES ? ######


poq <- mc.quant(rn.post, p = c(.9,.95,.99), lh = "gev")
#prq <- mc.quant(rn.prior, p = c(.1,.01,.001))
plot(density(poq[,1], adj = 2))
#lines(density(prq[,1], adj = 2), lty = 2)
plot(density(poq[,2], adj = 2))
#lines(density(prq[,2], adj = 2), lty = 2)
plot(density(poq[,3], adj = 2))
#lines(density(prq[,3], adj = 2), lty = 2)

# return levels
rl.pst(rn.post, lh = "gev")



##### Predictive distribution  ########

rl.pred(rn.prior, period = c(1,2,10), qlim = c(25,40), lh = "gev")
rl.pred(rn.post, period = c(1,2,10), qlim = c(30,40), lh = "gev")

###### Diagnostics : sensitivity analysis

reprn <- cbind(matrix(0, nrow = 116, ncol = 3), max_years$data)

for(i in 1:3) {
  j <- 2000*(i-1) + 1
  reprn[,i] <- evd::rgev(116, rn.post[j,1], rn.post[j,2], rn.post[j,3])
}
range(reprn) ; par(mfrow = c(2,2))
for(i in 1:4) hist(reprn[,i], freq = FALSE, breaks = seq(30,37,10))

# More formal way :

reprn <- matrix(0, nrow = 54, ncol = 8001)
for(i in 1:8001) {
  reprn[,i] <- evd::rgev(54, rn.post[i,1], rn.post[i,2], rn.post[i,3])  }
repmax <- apply(reprn, 2, max)
hist(repmax, freq = FALSE) ; abline(v=max(max_years$data), lwd = 3)
pv <- round(sum(repmax > max(max_years$data))/8001, 2)
text(300,.006, paste("p =", pv))
# The basic method of sensitivity analysis is to fit several models to the same problem.
# Posterior inferences from each model can then be compared. Posterior inferences
# will typically include marginal posterior distributions of the 3 parameters
# posterior distributions of GEV quantiles and posterior predictive distributions.
# The sensitivity of the marginal posterior density of the shape
# parameter ξ is often of particular interest.


######## Linear Trend , Variable threshold ########

# Scale and center the index t
# change values
shape <- c(38.9,7.1,47) ; scale <- c(1.5,6.3,2.6)
pr.quant <- prior.quant(shape = shape, scale = scale, trendsd = 5)
n <- 10000 ;   t0 <- c(32, 1.5,0,0) ;   s <- c(2,.35,.07,25) ;  b <- 2000
rn.prior2 <- posterior(n, t0, pr.quant, lh = "none", psd = s, burn = b)
t0 <- c(35,2,0,1) ;   s <- c(2,.2,.07,4) ; tt <- (1:20820 - 6576)/14610
rn.post2 <- posterior(n, t0, pr.quant, lh = "pp", data = max_all, thresh = 32,
                      noy = 54, trend = tt, psd = s, burn = b)


# Vary the threshold : introduce varying threshold

rl.pred()




# See Simulation-eficient shortest proba intervals PDF
library(quadprog)
library(SPIn)   # See source func
library(revdbayes)

vignette("revdbayes-vignette", package = "revdbayes")
# \url{https://cran.r-project.org/web/packages/revdbayes/vignettes/revdbayes-vignette.html}
# The essential difference between these two packages is that evdbayes performs sampling
# using Markov Chain Monte Carlo (MCMC) techniques, whereas revdbayes uses
# the generalised ratio-of-uniforms method (Wakefield, Gelfand, and Smith 1991),
# implemented using the rust package (Northrop 2016).
# https://cran.r-project.org/web/packages/rstan/vignettes/rstan.html