The fitdistcp package contains functions that generate predictive distributions
for various statistical models,
with and without parameter uncertainty.
Parameter uncertainty is included by using Bayesian prediction with a type of objective
prior known as a calibrating prior.
Calibrating priors are chosen to give predictions that give good reliability
(i.e., are well calibrated), for any
underlying true parameter values.
There are five functions for each model,
each of which uses training data x.
For model **** the five functions are as follows:
q****_cp returns predictive quantiles at the specified probabilities p,
and various other diagnostics.
r****_cp returns n random deviates from the predictive distribution.
d****_cp returns the predictive density function at the specified values y
p****_cp returns the predictive distribution function at the specified values y
t****_cp returns n random deviates from the posterior distribution
of the model parameters.
The q, r, d, p routines return two sets of results,
one based on maximum likelihood, and the other based on a calibrating prior.
The prior used depends on the model, and
is given under Details below.
Where possible, the Bayesian prediction integral is solved analytically. Otherwise, DMGS asymptotic expansions are used. Optionally, a third set of results is returned that integrates the prediction integral by sampling the parameter posterior distribution using the RUST rejection sampling algorithm.
qgumbel_cp(
x,
p = seq(0.1, 0.9, 0.1),
d1 = 0.01,
fd2 = 0.01,
means = FALSE,
waicscores = FALSE,
logscores = FALSE,
dmgs = TRUE,
rust = FALSE,
nrust = 1e+05,
debug = FALSE,
aderivs = TRUE
)rgumbel_cp(
n,
x,
d1 = 0.01,
fd2 = 0.01,
rust = FALSE,
mlcp = TRUE,
debug = FALSE,
aderivs = TRUE
)
dgumbel_cp(
x,
y = x,
d1 = 0.01,
fd2 = 0.01,
rust = FALSE,
nrust = 1000,
debug = FALSE,
aderivs = TRUE
)
pgumbel_cp(
x,
y = x,
d1 = 0.01,
fd2 = 0.01,
rust = FALSE,
nrust = 1000,
debug = FALSE,
aderivs = TRUE
)
tgumbel_cp(n, x, d1 = 0.01, fd2 = 0.01, debug = FALSE)
q**** returns a list containing at least the following:
ml_params: maximum likelihood estimates for the parameters.
ml_value: the value of the log-likelihood at the maximum.
standard_errors: estimates of the standard errors on the parameters,
from the inverse observed information matrix.
ml_quantiles: quantiles calculated using maximum likelihood.
cp_quantiles: predictive quantiles calculated using a calibrating prior.
maic: the AIC score for the maximum likelihood model, times -1/2.
cp_method: a comment about the method used to generate the cp prediction.
For models with predictors, q**** additionally returns:
predictedparameter: the estimated value for parameter,
as a function of the predictor.
adjustedx: the detrended values of x
r**** returns a list containing the following:
ml_params: maximum likelihood estimates for the parameters.
ml_deviates: random deviates calculated using maximum likelihood.
cp_deviates: predictive random deviates calculated using a calibrating prior.
cp_method: a comment about the method used to generate the cp prediction.
d**** returns a list containing the following:
ml_params: maximum likelihood estimates for the parameters.
ml_pdf: density function from maximum likelihood.
cp_pdf: predictive density function calculated using a calibrating prior
(not available in EVT routines, for mathematical reasons, unless using RUST).
cp_method: a comment about the method used to generate the cp prediction.
p*** returns a list containing the following:
ml_params: maximum likelihood estimates for the parameters.
ml_cdf: distribution function from maximum likelihood.
cp_cdf: predictive distribution function calculated using a calibrating prior
(not available in EVT routines, for mathematical reasons, unless using RUST).
cp_method: a comment about the method used to generate the cp prediction.
t*** returns a list containing the following:
theta_samples: random samples from the parameter posterior.
a vector of training data values
a vector of probabilities at which to generate predictive quantiles
if aderivs=FALSE, the delta used for numerical derivatives
with respect to the first parameter
if aderivs=FALSE, the fractional delta used for numerical derivatives
with respect to the second parameter
logical that indicates whether to run additional calculations and return analytical estimates for the distribution means (longer runtime)
logical that indicates whether to run additional calculations and return estimates for the WAIC1 and WAIC2 scores (longer runtime)
logical that indicates whether to run additional calculations and return leave-one-out estimates of the log-score (much longer runtime, non-EVT models only)
logical that indicates whether DMGS calculations should be run or not (longer run time)
logical that indicates whether RUST-based posterior sampling calculations should be run or not (longer run time)
the number of posterior samples used in the RUST calculations
logical for turning on debug messages
(for code testing only) logical for whether to use analytic derivatives (instead of numerical). By default almost all models now use analytical derivatives.
the number of random samples required
logical that indicates whether maxlik and parameter uncertainty calculations should be performed (turn off to speed up RUST)
a vector of values at which to calculate the density and distribution functions
The Gumbel distribution has distribution function $$F(x;\mu,\sigma)=\exp\left(-\exp\left(-\frac{x-\mu}{\sigma}\right)\right)$$ where \(x\) is the random variable and \(\mu,\sigma>0\) are the parameters.
The calibrating prior is given by the right Haar prior, which is $$\pi(\sigma) \propto \frac{1}{\sigma}$$ as given in Jewson et al. (2025).
q**** optionally returns the following:
If rust=TRUE:
ru_quantiles: predictive quantiles calculated using a calibrating prior,
using posterior sampling with the RUST algorithm, based on inverting
an empirical CDF based on nrust samples.
If waicscores=TRUE:
waic1: the WAIC1 score for the calibrating prior model.
waic2: the WAIC2 score for the calibrating prior model.
If logscores=TRUE:
ml_oos_logscore: the leave-one-out logscore for the
maximum likelihood prediction
(not available in EVT routines, for mathematical reasons)
cp_oos_logscore: the leave-one-out logscore for the
parameter uncertainty model
available in EVT routines, for mathematical reasons)
If means=TRUE:
ml_mean: analytic estimate of the mean of the MLE predictive distribution,
where possible
cp_mean: analytic estimate of the mean of the calibrating prior
predictive distribution, where mathematically possible.
Can be compared with the mean estimated from random deviates.
r**** optionally returns the following:
If rust=TRUE:
ru_deviates: nrust predictive random deviatives calculated using a
calibrating prior,
using posterior sampling with RUST.
d**** optionally returns the following:
If rust=TRUE:
ru_pdf: predictive density calculated using a calibrating prior,
using posterior sampling with RUST,
averaging over nrust density functions.
p**** optionally returns the following:
If rust=TRUE:
ru_cdf: predictive probability calculated using a calibrating prior,
using posterior sampling with RUST,
averaging over nrust distribution functions.
Selecting these additional outputs increases runtime. They are optional so that runtime for the basic outputs is minimised. This facilitates repeated experiments that evaluate reliability over many thousands of repeats.
This model is a homogeneous model, and the cp results are based on the
right Haar prior.
For homogeneous models
(models with sharply transitive transformation groups),
a Bayesian prediction based on the right Haar prior
gives exact reliability,
as shown by Severini et al. (2002),
even when the true parameters are unknown.
This means that probabilities in the prediction will correspond to frequencies
of future outcomes in repeated trials (if the model is correct).
Maximum likelihood prediction does not give reliable predictions, even when the model is correct, because it does not account for parameter uncertainty. In particular, maximum likelihood predictions typically underestimate the tail in repeated trials.
The reliability of the maximum likelihood and the calibrating prior
predictive quantiles produced by
the q****_cp routines in fitdistcp can be quantified
using repeated simulations with the routine reltest.
For this model, the Bayesian prediction equation cannot be solved analytically,
and is approximated using the DMGS asymptotic
expansions given by Datta et al. (2000).
This approximation seems to work well for medium and large sample sizes,
but may not work well for small sample sizes (e.g., <20 data points).
For small sample sizes, it may be preferable to use posterior
sampling by setting rust=TRUE and looking at the ru outputs.
The performance for any sample size, in terms of reliability,
can be tested using reltest.
The Bayesian prediction equation can also be integrated using
ratio-of-uniforms-sampling-with-transformation (RUST), using the option
rust=TRUE.
fitdistcp then calls Paul Northrop's rust package (Northrop, 2023).
The RUST calculations are slower than the DMGS calculations.
For small sample sizes (e.g., n<20), and the very extreme tail, the DMGS approximation is somewhat poor (although always better than maximum likelihood) and it may be better to use RUST. For medium sample sizes (30+), DMGS is reasonably accurate, except for the very far tail.
It is advisable to check the RUST results for convergence versus the number of RUST samples.
It may be interesting to compare the DMGS and RUST results.
Stephen Jewson stephen.jewson@gmail.com
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,
#
# example 1
x=fitdistcp::d50gumbel_example_data_v1
p=c(1:9)/10
q=qgumbel_cp(x,p,rust=TRUE,nrust=1000)
xmin=min(q$ml_quantiles,q$cp_quantiles);
xmax=max(q$ml_quantiles,q$cp_quantiles);
plot(q$ml_quantiles,p,xlab="quantile estimates",xlim=c(xmin,xmax),
sub="(from qgumbel_cp)",
main="Gumbel: quantile estimates");
points(q$cp_quantiles,p,col="red",lwd=2)
points(q$ru_quantiles,p,col="blue")
Run the code above in your browser using DataLab