Fitdistcp: fitting statistical models using calibrating priors

fitdistcp is a free R package for fitting statistical models using calibrating priors, with the goal of making reliable predictions.

fitdistcp can be called from R or Python (details below), and we hope to release a pure Python version in summer 2025.

Maximum likelihood and other point estimate methods underestimate predictive tail probabilities. For instance, with maximum likelihood prediction, the GEVD with one predictor and shape parameter of -0.3 fitted to 50 data points underestimates the 1 in 50 probability by a factor of 4. This is called a reliability bias.

Calibrating prior prediction provides an objective way to reduce this reliability bias by including parameter uncertainty. It eliminate the bias completely in many cases. The benefit of calibrating prior prediction, relative to maximum likelihood, is largest for small samples, models with more parameters, and in the tails.

fitdistcp uses objective Bayesian methods. The prior is chosen to maximize the reliability of the predictions, i.e., to give predictive probabilities that are as close as possible to the frequencies of future outcomes. The integration is performed analytically if possible, or using DMGS asymptotics otherwise. Posterior sampling is also available as an option (but is slower).

For many cases, charts are provided below that quantify the size of the maximum likelihood reliability bias, and the reduced bias from calibrating prior prediction integrated with DMGS.

Models

This is a list of the models that are currently supported, in alphabetical order. This list is motivated by various climate and actuarial applications. Let me know if you have any suggestions for other models to include.

• The R command needs to be prefixed with one of d,p,q,r,t
• H/I indicates whether the model is homogenenous or inhomogeneous
• A/D indicates whether the default integration method is analytic or DMGS
• In some cases, reliability bias charts are provided for maximum likelihood and calibrating prior predictions.
• For EVT models, the reliability bias charts are given as a function of the true shape parameter

Examples

• Install fitdistcp from CRAN using install.packages("fitdistcp").
• Load fitdistcp using library(fitdistcp).

Example 1: Fitting a normal distribution and calculating the fitted probability density function

> x=rnorm(20)		# make some example training data 
> y=seq(-5,5,0.01)	# define the points at which we wish to calculate the density
> d=dnorm_cp(x,y)	# this command calculates two sets of predictive densities for the normal, one based on maxlik, 
			# and one that includes parameter uncertainty based on a calibrating prior
> print(d$ml_pdf)	# have a look at the maxlik densities
> print(d$cp_pdf)	# have a look at the densities that include parameter uncertainty, which will have fatter tails

					

Example 2: Fitting a GEV distribution and calculating the fitted quantiles

> set.seed(1)
> x=rlogis(20)				# make some example training data 
> p=seq(0.001,0.999,0.001)		# define the probabilities at which we wish to calculate the quantiles
> q=qgev_cp(x,p)			# this command calculates two sets of predictive quantiles for the GEV, 
					# one based on maxlik, and one that includes parameter uncertainty based on a calibrating prior
> q$ml_params 				# have a look at the maxlik parameters
> plot(q$ml_quantiles,p)		# plot the maxlik quantiles
> points(q$cp_quantiles,p,col="red")	# overplot the quantiles that include parameter uncertainty, which will have fatter tails

					

Example 3: Fitting a GPD and calculating the fitted cumulative distribution function

> set.seed(1)
> x=rlnorm(20,sdlog=2)		# make some example training data (all positive values)
> y=seq(0.1,5,0.1)		# define the points at which we wish to calculate the probabilities
> p=pgpd_k1_cp(x,y)		# this command calculates two sets of predictive probabilities for the GPD (with known location parameter, hence k1), 
				# one based on maxlik, and one that includes parameter uncertainty based on a calibrating prior
> plot(y,p$ml_cdf)		# have a look at the maxlik probabilities
> point(y,p$cp_pdf,col="red")	# have a look at the probabilities that include parameter uncertainty, 
				# which will have fatter tails

					

Example 4: Fitting a Gumbel distribution and generating random numbers

> x=rnorm(10)					# make some example training data 
> r=rgumbel_cp(10000,x)				# this command calculates two sets of 10,000 predictive random numbers, 
						# based on maxlik, and then with parameter uncertainty based on a calibrating prior
> plot(density(r$ml_deviates))			# plot the density of the maxlik random numbers
> lines(density(r$cp_deviates),col="red")	# overplot the density of the random numbers including parameter uncertainty,
						# which will have fatter tails

					

Example 5: Generating posterior samples from a Weibull distribution

> x=rlnorm(20)						# make some example training data (all positive values)
> t=tweibull_cp(1000,x)$theta_samples			# this command generates 1000 posterior samples
			

Example 6: Fitting a GEV distribution with a predictor on the location parameter

> set.seed(1)
> n=30
> t=0.1*c(1:n)/n			# make an example predictor
> x=t+rlogis(n)				# make training data based on the predictor
> t0=t[n]				# set the value of the predictor at which to make the prediction
> p=seq(0.001,0.999,0.001)		# specify the probabilities at which to generate quantiles
> q=qgev_p1_cp(x,t,t0=t0,n0=NA,p) 	# specify either t0, or n0, and make the other one NA
			

Example 7: Fitting a GEV distribution with predictors on location and scale parameters

> set.seed(1)
> n=50
> t1=0.1*c(1:n)/n					# make the first predictor
> t2=0.2*c(1:n)/n					# make the second predictor
> x=(t1+rlogis(n))*t2					# make training data based on the predictors
> n0=n							# set the values of the predictors at which to make the prediction
> p=seq(0.001,0.999,0.001)		
> q=qgev_p12_cp(x,t1,t2,t01=NA,t02=NA,n01=n0,n02=n0,p)	# specify either t0, or n0, and make the other one NA
			

Detailed Examples by model

Theory

• The goal is to make predictions that are as reliable as possible. Reliable means that events that are predicted to occur with a probability of 1% will indeed occur 1% of the time. Predictions from point estimate methods such as maximum likelihood are not reliable, and underestimate the tails.
• The predictions with parameter uncertainty are generated using the Bayesian prediction equation.
• Bayesian prediction requires a prior. For homogeneous models (i.e., models with sharply transitive transformation groups), we use the right Haar prior, which gives predictions that are in principle exactly reliable (if the model assumptions are correct, and depending on the integration scheme used). This applies to our Cauchy, exponential, Frechet, generalized normal, GEV with known shape, Gumbel, halfnormal, logistic, log-normal, normal, Pareto, t and Weibull models. It also applies to all these models with predictors.
• For inhomogeneous models (i.e., models without sharply transitive transformation groups) there is no solution that gives exact reliability, so we use simple priors that we have found will give close to exact reliability. They are in all cases more reliable than maxlik, but may not be very reliable for small sample size (e.g., n<50). This applies to GEV, GPD, gamma, inverse gamma and inverse Gaussian, and these models with predictors.
• Simulations can be used to test whether a certain prior gives reliable predictions, as a function of the true parameters. Simulation results are given in the reference paper cited below, and can be reproduced using the command reltest.
• Where possible, we integrate the Bayesian prediction equation analytically. This applies to the exponential, normal, log-normal, uniform and Pareto.
• Otherwise, by default we integrate the Bayesian prediction equation using the DMGS asymptotic expansion, which works well except for very small samples (n<10) (although even for very small sample sizes calibrating priors with DMGS still give more reliable predictions than maximum likelihood). DMGS is very fast, especially for quantiles (~1000x faster than MCMC).
• DMGS requires the calculation of various derivatives. Where possible, we use analytic expressions. Otherwise, we use simple numerical derivatives.
• There is also an option to integrate using posterior sampling, which then uses Paul Northrop's rust library (ratio of uniforms with transformation). Posterior sampling is slower than DMGS for the same level of accuracy (and for quantiles, much slower), but is ultimately more accurate if enough posterior samples are used.

Other

• All calculations return the maxlik parameter estimates
• All calculations with the q***_cp routines calculate the AIC and standard errors. AIC is presented as MAIC, which is -0.5 times the standard AIC, for better comparison with log scores, and so that higher is better.
• Optionally, WAIC and out-of-sample log-score can be calculated too (although oos log-score can't be calculated for EVT routines, for mathematical reasons)
• Model selection can be illustrated using the ms_flat_1tail, ms_flat_2tail, ms_predictors_1tail and ms_predictors_2tail routines.

Python

> import numpy as np 					# regular Python imports
> import os			
> os.environ["R_HOME"] = "C:/PROGRA~1/R/R-43~1.1"	# point to your R installation
> import rpy2.robjects as robjects			# import the routines that make R work in Python
> from rpy2.robjects.packages import importr
> cp=importr('fitdistcp')				# define the R library we want to use
> x = np.arange(1, 11) ** 4				# define the input data, as Python variables
> p = np.arange(1, 4) / 4				# noting that Python counts from 0 not 1
> rx=robjects.FloatVector(x)				# convert the input data to r variables
> rp=robjects.FloatVector(p)
> q=cp.qexp_cp(rx,rp)					# call the R function by prefixing with the library name we defined above
> ml=q.rx2("ml_quantiles")				# extract the parts of the output that we want, into Python variables
> cp=q.rx2("cp_quantiles")
> print("ml quantiles:")				# have a look at the results
> print(ml)
> print("cp quantiles:")
> print(cp)
> # robjects.r.source("script.R", encoding="utf-8")	# how to run the whole thing as a script
					

Citation

• S. Jewson, T. Sweeting and L. Jewson (2024): Reducing Reliability Bias in Assessments of Extreme Weather Risk using Calibrating Priors: ASCMO (Advances in Statistical Climatology, Meteorology and Oceanography)
• Bibtex:

  @article{jewsonet2024,
  	author = {Jewson, Stephen and Sweeting, Trevor and Jewson, Lynne},
  	title = {Reducing reliability bias in assessments of extreme weather risk using calibrating priors},
  	journal = {ASCMO},
  	volume = {11},
  	number = {1},
  	pages = {1-22},
  	doi = {10.5194/ascmo-11-1-2025},
  	url = {https://doi.org/10.5194/ascmo-11-1-2025},
  	year = {2025}
  }
  							

Notes