fitdistcp PackageUses simulations to evaluate the reliability of
the predictive quantiles produced by the q****_cp routines in the fitdistcp package.
reltest(
model = "exp",
ntrials = 1000,
nrepeats = 3,
nx = 20,
params = c(1),
alpha = seq(0.005, 0.995, 0.005),
plotflag = TRUE,
verbose = TRUE,
dmgs = TRUE,
debug = FALSE,
aderivs = TRUE,
unbiasedv = FALSE,
pwm = FALSE,
minxi = -10,
maxxi = 10
)A plot showing 9 different reliability checks, and a list containing various outputs, including the probabilities shown in the plot.
which distribution to test. Possibles values are
"exp",
"pareto_k1",
"halfnorm",
"unif",
"norm",
"norm_dmgs",
"gnorm_k3",
"lnorm",
"lnorm_dmgs",
"logis",
"lst_k3",
"cauchy",
"gumbel",
"frechet_k1",
"weibull",
"gev_k3",
"exp_p1",
"pareto_p1k3",
"norm_p1",
"lnorm_p1",
"logis_p1",
"lst_p1k4",
"cauchy_p1",
"gumbel_p1",
"frechet_p2k1",
"weibull_p2",
"gev_p1k4",
"norm_p12",
"lst_p12k5",
"gamma",
"invgamma",
"invgauss",
"gev",
"gpd_k1",
"gev_p1".
"gev_p12".
"gev_p123".
the number of trials to run. 5000 typically gives good results.
the number of entire repeats of the test to run, to check for convergence. 3 is a good choice.
the length of the training data to use.
values for the parameters for the specified distribution
the exceedance probability values at which to test
logical to turn the plotting on and off
logical to turn loop counting on and off
logical to turn DMGS calculations on and off (to optimize speed for maxlik only calculations)
logical for turning debug messages on and off
logical for whether to use analytic derivatives (instead of numerical)
logical for whether to use the unbiased variance instead of maxlik (for the normal)
logical for whether to use PWM instead of maxlik (for the GEV)
minimum value for EVT shape parameter
maximum value for EVT shape parameter
Stephen Jewson stephen.jewson@gmail.com
The maximum likelihood quantiles (plotted in blue) do not give good reliability. They typically underestimate the tails (see panel (f)).
For
"exp",
"pareto_k1",
"unif",
"norm",
"lnorm",
"norm_p1" and
"lnorm_p1",
the calibrating prior quantiles are calculated using the right Haar prior
and an exact solution for the Bayesian prediction integral.
They will converge towards exact reliability with a large enough number of trials,
for any sample size.
For
"halfnorm",
"norm_dmgs",
"lnorm_dmgs",
"gnorm_k3",
"logis",
"lst_k3",
"cauchy",
"gumbel",
"frechet_k1",
"weibull",
"gev_k3",
"exp_p1",
"pareto_p1k3",
"gumbel_p1",
"logis_p1" and
"lst_p1k4"
"cauchy_p1",
"gumbel_p1",
"frechet_p2k1",
"weibull_p2",
"gev_p1k4",
"norm_p12",
"lst_p12k5"
the calibrating prior quantiles are calculated using the right Haar prior,
with the DMGS asymptotic solution for the Bayesian prediction integral.
They will converge towards good reliability with a large enough number of trials,
with the only deviation from exact reliability being due to the neglect of
higher order terms in the asymptotic expansion.
They will converge towards exact reliability with a large enough number of trials
and a large enough sample size.
For
"gamma",
"invgamma",
"invgauss",
"gev",
"gpd_k1" and
"gev_p1",
"gev_p12",
"gev_p123",
the calibrating prior quantiles are calculated using the "fitdistcp"
recommended calibrating priors,
with the DMGS asymptotic solution for the Bayesian prediction integral.
The chosen priors give reasonably good reliability with a
large enough number of trials,
and for large sample sizes, but may give poor reliability for small
sample sizes (e.g., n<20).
If you use this package, we would be grateful if you would cite the following reference, which gives the various calibrating priors, and tests them for reliability:
Jewson S., Sweeting T. and Jewson L. (2024): Reducing Reliability Bias in Assessments of Extreme Weather Risk using Calibrating Priors; ASCMO Advances in Statistical Climatology, Meteorology and Oceanography), https://ascmo.copernicus.org/articles/11/1/2025/.
An introduction to fitdistcp, with more examples,
is given on this webpage.
The fitdistcp package currently includes the following models (in alphabetical order):
Cauchy (cauchy),
Cauchy with linear predictor on the mean (cauchy_p1),
Exponential (exp),
Exponential with log-linear predictor on the scale (exp_p1),
Frechet with known location parameter (frechet_k1),
Frechet with log-linear predictor on the scale and known location parameter
(frechet_p2k1),
Gamma (gamma),
Generalized normal (gnorm),
GEV (gev),
GEV with linear predictor on the location (gev_p1),
GEV with linear predictor on the location and log-linear prediction on the scale
(gev_p12),
GEV with linear predictor on the location, log-linear prediction on the scale,
and linear predictor on the shape (gev_p123),
GEV with linear predictor on the location and known shape (gev_p1k3),
GEV with known shape (gev_k3),
GPD with known location (gpd_k1),
Gumbel (gumbel),
Gumbel with linear predictor on the mean(gumbel_p1),
Half-normal (halfnorm),
Inverse gamma (invgamma),
Inverse Gaussian (invgauss),
t distribution with unknown location and scale and known DoF (lst_k3),
t distribution with unknown location and scale, linear predictor on the location,
and known DoF (lst_p1k3),
Logistic (logis),
Logistic with linear predictor on the location (logis_p1),
Log-normal (lnorm),
Log-normal with linear predictor on the location (lnorm_p1),
Normal (norm),
Normal with linear predictor on the mean (norm_p1),
Pareto with known scale (pareto_k2),
Pareto with log-linear predictor on the shape and known scale (pareto_p1k2),
Uniform (unif),
Weibull (weibull),
Weibull with linear predictor on the scale (weibull_p2),
The level of predictive probability matching achieved
by the maximum likelihood and calibrating prior quantiles, for any model,
sample size and true parameter values, can be demonstrated using the
routine reltest.
Model selection among models can be demonstrated using the routines
ms_flat_1tail,
ms_flat_2tail,
ms_predictors_1tail,
and ms_predictors_2tail,
set.seed(1)
# example 1
# -runs the default settings, which test reliability for the exponential distribution
reltest()
Run the code above in your browser using DataLab