W.diag | R Documentation |
Function to produce diagnostic plots and test statistics for the threshold diagnostics exploiting structure of maximum likelihood estimators based on the non-homogeneous Poisson process likelihood
W.diag(
xdat,
model = c("nhpp", "exp", "invexp"),
u = NULL,
k,
q1 = 0,
q2 = 1,
par = NULL,
M = NULL,
nbs = 1000,
alpha = 0.05,
plots = c("LRT", "WN", "PS"),
UseQuantiles = FALSE,
changepar = TRUE,
...
)
xdat |
a numeric vector of data to be fitted. |
model |
string specifying whether the univariate or bivariate diagnostic should be used. Either |
u |
optional; vector of candidate thresholds. |
k |
number of thresholds to consider (if |
q1 |
lowest quantile for the threshold sequence. |
q2 |
upper quantile limit for the threshold sequence ( |
par |
parameters of the NHPP likelihood. If |
M |
number of superpositions or 'blocks' / 'years' the process corresponds to (can affect the optimization) |
nbs |
number of simulations used to assess the null distribution of the LRT, and produce the p-value |
alpha |
significance level of the LRT |
plots |
vector of strings indicating which plots to produce; |
UseQuantiles |
logical; use quantiles as the thresholds in the plot? |
changepar |
logical; if |
... |
additional parameters passed to |
The function is a wrapper for the univariate (non-homogeneous Poisson process model) and bivariate exponential dependence model.
For the latter, the user can select either the rate or inverse rate parameter (the inverse rate parametrization works better for uniformity
of the p-value distribution under the LR
test.
There are two options for the bivariate diagnostic: either provide pairwise minimum of marginally
exponentially distributed margins or provide a n
times 2 matrix with the original data, which
is transformed to exponential margins using the empirical distribution function.
plots of the requested diagnostics and an invisible list with components
MLE
: maximum likelihood estimates from all thresholds
Cov
: joint asymptotic covariance matrix for \xi
, \eta
or 1/\eta
.
WN
: values of the white noise process
LRT
: values of the likelihood ratio test statistic vs threshold
pval
: P-value of the likelihood ratio test
k
: final number of thresholds used
thresh
: threshold selected by the likelihood ratio procedure
qthresh
: quantile level of threshold selected by the likelihood ratio procedure
cthresh
: vector of candidate thresholds
qcthresh
: quantile level of candidate thresholds
mle.u
: maximum likelihood estimates for the selected threshold
model
: model fitted
Jennifer L. Wadsworth
Wadsworth, J.L. (2016). Exploiting Structure of Maximum Likelihood Estimators for Extreme Value Threshold Selection, Technometrics, 58(1), 116-126, http://dx.doi.org/10.1080/00401706.2014.998345
.
## Not run:
set.seed(123)
# Parameter stability only
W.diag(xdat = abs(rnorm(5000)), model = 'nhpp',
k = 30, q1 = 0, plots = "PS")
W.diag(rexp(1000), model = 'nhpp', k = 20, q1 = 0)
xbvn <- mvrnorm(n = 6000,
mu = rep(0, 2),
Sigma = cbind(c(1, 0.7), c(0.7, 1)))
# Transform margins to exponential manually
xbvn.exp <- -log(1 - pnorm(xbvn))
#rate parametrization
W.diag(xdat = apply(xbvn.exp, 1, min), model = 'exp',
k = 30, q1 = 0)
W.diag(xdat = xbvn, model = 'exp', k = 30, q1 = 0)
#inverse rate parametrization
W.diag(xdat = apply(xbvn.exp, 1, min), model = 'invexp',
k = 30, q1 = 0)
## End(Not run)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.