Create a Posterior for a Poisson-GP Model

Create an object representing a "Poor Man's" Posterior for a Poisson-GP model, mainly using MCMC iterates.

Usage

poisGPBayes0(MCMC, threshold,
             data = NULL, effDuration,
             MAX.data = NULL, MAX.effDuration = NULL,
             OTS.data = NULL, OTS.threshold = NULL, OTS.effDuration = NULL,
             MAP = NULL,
             nOT = NA,
             a0 = 1.0, b0 = 0.0)

Arguments

MCMC: An object that can be coerced into a matrix containing the MCMC iterates. It should have the burnin period removed and be thinned if necessary.
threshold: the POT threshold.
data: An optional structure containing the observations Over the Threshold. The user should make sure that these data are identical to those used to obtain MCMC. It can be a numeric vector or an object inheriting from one of the two classes "Rendata" (from Renext) and "potData" (from potomax) in which case it will be coerced into the class "potData".
effDuration: The observation effective duration. When data can be coerced into potData, the provided value for effDuration will be ignored and the value of object$OT$effDuration will be used instead.
MAX.data, MAX.effDuration: Optional censored data (type block maxima or $r$-largest) as used in poisGP. Ignored if data inherits from one of the two classes "Rendata" and "potData"
OTS.data, OTS.threshold, OTS.effDuration: Optional censored data as used in poisGP. Ignored if data inherits from one of the two classes "Rendata" and "potData"
MAP: An optional vector of Maximum A Posteriori for the parameter vector. Should be named with names matching the colnames of MCMC.
nOT: An optional integer giving the number of exceedances during the observation period. This number is needed when MCMC does not have a column for the rate of the Poisson process. It should be consistent with MCMC.
a0, b0: The shape and the rate of the prior for lambda (used only when MCMC has two columns).

Value

An object with class "GEVBayes0" inheriting from "Bayes0". This object can be used to produce RL plots.

Details

The Poisson-GP model involves three parameters: one is the rate lambda of the Poisson process describing the exceedances over the threshold. The two other parameters are the scale and the shape of the GP Distribution used for the excesses over the threshold. The object MCMC can have either three columns or two columns only. In the second case, the column corresponding to the rate lambda is omitted; it is then assumed that lambda is a posteriori independent of the GP parameters and that it has a gamma posterior margin. The shape of the posterior distribution is $a_n := a_0 + n_{\textrm{OT}}$ and the rate is $b_n := b_0 + w_{\textrm{OT}}$, where $w_{\textrm{OT}}$ and $w_{\textrm{OT}}$ are the number of exceedances and the duration of the observation period.

Caution

The user must be sure that the data provided by using the dedicated arguments are consistent with the MCMC iterates provided in MCMC. The 'data' dedicated arguments are: data, effDuration or/and MAX.data, MAX.effDuration OTS.data, OTS.threshold, OTS.effDuration. Since attaching data to the object is essentially for graphical purpose, it can be simpler not to attach the data and to create a separated object with class "potData". This object can be used with its autolayer method.

Examples

require(revdbayes)
## ========================================================================
## use the 'rainfall data': daily rainfalls over 57 years
## ========================================================================
data(rainfall)
rainfall2 <- rainfall[!is.na(rainfall)]
u <- 40  ## threshold (mm)
w <- 57  ## obs duration (year)
nOT <- sum(rainfall2 > u)
prior <- set_prior(prior = "flatflat", model = "gp")
nSim <- 10000

## ========================================================================
## get a posterior for the "excess part" of the model
## ========================================================================
postGP <- rpost(n = 10000, model = "gp", prior = prior,
                data = rainfall2, thresh = u)
MCMCGP <- postGP$sim_vals
MAP <- postGP$f_mode
names(MAP) <- c("scale", "shape")

## ========================================================================
## assuming a "flat gamma prior" for 'lambda', add MCMC iterates
## for 'lambda' to those existing for GP (new matrix column)
## ========================================================================
MCMCpoisGP <- cbind(lambda = rgamma(nSim, shape = 1 + nOT, rate = 0 + w),
                    MCMCGP)

## ========================================================================
## create the object
## ========================================================================
post0 <- poisGPBayes0(MCMC = MCMCpoisGP, threshold = u, effDuration = w) 
summary(post0)
#> poisGP Model Bayesian Inference
#> o Number of OT observations: NA
#> o Number of MCMC iterates: 10000
#> o Posterior mean [sd]:
#>      lambda       scale       shape 
#> 1.53 [0.16] 7.99 [1.27] 0.19 [0.13] 
coef(post0)
#>          lambda    scale     shape
#> mean   1.526300 7.987665 0.1947722
#> median 1.522316 7.907451 0.1802754
#> mode         NA       NA        NA

## ========================================================================
## 'lambda' and the GP parameters are a posteriori independent
## ========================================================================
cov2cor(vcov(post0))
#>              lambda        scale        shape
#> lambda  1.000000000 -0.003300372  0.003631015
#> scale  -0.003300372  1.000000000 -0.593252985
#> shape   0.003631015 -0.593252985  1.000000000

## ========================================================================
## create an object assuming 'lambda' independent of the GP param
## ========================================================================
post1 <- poisGPBayes0(MCMC = MCMCGP, threshold = u,
                      effDuration = w, nOT = nOT, MAP = MAP)