<figure><img src="image_1.png" /></figure>

HAL
open science

# fitdistrplus: An R Package for Fitting Distributions

Marie Laure Delignette-Muller, Christophe Dutang

► To cite this version:

Marie Laure Delignette-Muller, Christophe Dutang. fitdistrplus: An R Package for Fitting Distributions. Journal of Statistical Software, 2015, 64 (4), 10.18637/jss.v064.i04. hal-01616147

HAL Id: hal-01616147

https://hal.science/hal-01616147v1

Submitted on 11 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# fitdistrplus: An R Package for Fitting Distributions

Marie Laure Delignette-Muller

Université de Lyon

Christophe Dutang

Université de Strasbourg

## Abstract

The package **fitdistrplus** provides functions for fitting univariate distributions to different types of data (continuous censored or non-censored data and discrete data) and allowing different estimation methods (maximum likelihood, moment matching, quantile matching and maximum goodness-of-fit estimation). Outputs of **fitdist** and **fitdistcens** functions are S3 objects, for which specific methods are provided, including **summary**, **plot** and **quantile**. This package also provides various functions to compare the fit of several distributions to the same data set and can handle to bootstrap parameter estimates. Detailed examples are given in food risk assessment, ecotoxicology and insurance contexts.

**Keywords:** probability distribution fitting, bootstrap, censored data, maximum likelihood, moment matching, quantile matching, maximum goodness-of-fit, distributions, R.

## 1. Introduction

Fitting distributions to data is a very common task in statistics and consists in choosing a probability distribution modeling the random variable, as well as finding parameter estimates for that distribution. This requires judgment and expertise and generally needs an iterative process of distribution choice, parameter estimation, and quality of fit assessment. In the R (R Core Team 2014) package MASS (Venables and Ripley 2010), maximum likelihood estimation is available via the **fitdistr** function; other steps of the fitting process can be done using other R functions (Ricci 2005). In this paper, we present the R package **fitdistrplus** (Delignette-Muller, Pouillot, Denis, and Dutang 2015) implementing several methods for fitting univariate parametric distributions. A first objective in developing this package was to provide R users with a set of functions dedicated to help this overall process.

The **fitdistr** function in the MASS package estimates distribution parameters by maximizing the likelihood function using the optim function. No distinction between parameters with different roles (e.g., main parameter and nuisance parameter) is made, as this paper focuses on parameter estimation from a general point-of-view. In some cases, other estima-
tion methods could be preferred, such as maximum goodness-of-fit estimation (also called minimum distance estimation), as proposed in the R package **actuar** with three different goodness-of-fit distances (Dutang, Goulet, and Pigeon 2008). While developing the **fitdistrplus** package, a second objective was to consider various estimation methods in addition to maximum likelihood estimation (MLE). Functions were developed to enable moment matching estimation (MME), quantile matching estimation (QME), and maximum goodness-of-fit estimation (MGE) using eight different distances. Moreover, the **fitdistrplus** package offers the possibility to specify a user-supplied function for optimization, useful in cases where classical optimization techniques, not included in **optim**, are more adequate.

In applied statistics, it is frequent to have to fit distributions to censored data (Klein and Moeschberger 2003; Helsel 2005; Busschaert, Geeraerd, Uyttendaele, and VanImpe 2010; Leha, Beissbarth, and Jung 2011; Commeau, Parent, Delignette-Muller, and Cornu 2012). The *fitdistr* function in the **MASS** package does not enable maximum likelihood estimation with this type of data. Some packages can be used to work with censored data, especially survival data (Therneau 2014; Hirano, Clayton, and Upper 1994; Jordan 2005), but those packages generally focus on specific models, enabling the fit of a restricted set of distributions. A third objective is thus to provide R users with a function to estimate univariate distribution parameters from right-, left- and interval-censored data.

Few packages on the Comprehensive R Archive Network (CRAN, http://CRAN.R-project.org) provide estimation procedures for any user-supplied parametric distribution and support different types of data. The **distrMod** package (Kohl and Ruckdeschel 2010) provides an object-oriented (S4) implementation of probability models and includes distribution fitting procedures for a given minimization criterion. This criterion is a user-supplied function which is sufficiently flexible to handle censored data, yet not in a trivial way, see Example M4 of the **distrMod** vignette. The fitting functions **MLEstimator** and **MDEstimator** return an S4 class for which a coercion method to class 'mle' is provided so that the respective functionalities (e.g., **confint** and **logLik**) from package **stats4** are available, too. In **fitdistrplus**, we chose to use the standard S3 class system because of its ease of understanding by most R users. When designing the **fitdistrplus** package, we did not forget to implement method functions for the S3 classes. Finally, various other packages provide functions to estimate the mode, the moments or the L-moments of a distribution, see the reference manuals of the packages **modeest** (Poncet 2012), **lmomco** (Asquith 2014) and **Lmoments** (Karvanen 2006, 2014).

This paper reviews the various features of **fitdistrplus** (Delignette-Muller et al. 2015). The package is available from the Comprehensive R Archive Network (CRAN) at http://CRAN.R-project.org/package=fitdistrplus. The development version of the package is located at R-Forge as one package of the project "Risk Assessment with R" (http://R-Forge.R-project.org/projects/riskassessment/). The paper is organized as follows: Section 2 presents tools for fitting continuous distributions to classic non-censored data. Section 3 deals with other estimation methods and other types of data, before Section 4 concludes.

## 2. Fitting distributions to continuous non-censored data

### 2.1. Choice of candidate distributions

For illustrating the use of various functions of the **fitdistrplus** package with continuous non-
<figure><img src="image_4.png" /><figcaption>Figure 1: Histogram and CDF plots of an empirical distribution for a continuous variable (serving size from the groundbeef data set) as provided by the plotdist function.</figcaption></figure>

censored data, we will first use a data set named groundbeef which is included in our package. This data set contains pointwise values of serving sizes in grams, collected in a French survey, for ground beef patties consumed by children under 5 years old. It was used in a quantitative risk assessment published by Delignette-Muller, Cornu, and AFSSA-STEC-Study-Group (2008).

```R
R> library("fitdistrplus")
R> data("groundbeef", package = "fitdistrplus")
R> str(groundbeef)
```
```
/data.frame:      254 obs. of 1 variable:
$ serving: num  30  10  20  24  20  24  40  20  50  30 ...
```

Before fitting one or more distributions to a data set, it is generally necessary to choose good candidates among a predefined set of distributions. This choice may be guided by the knowledge of stochastic processes governing the modeled variable, or, in the absence of knowledge regarding the underlying process, by the observation of its empirical distribution. To help the user in this choice, we developed functions to plot and characterize the empirical distribution.

First of all, it is common to start with plots of the empirical distribution function and the histogram (or density plot), which can be obtained with the plotdist function of the fitdistrplus package. This function provides two plots (see Figure 1): the left-hand plot is by default the histogram on a density scale (or density plot of both, according to values of arguments histo and demp) and the right-hand plot the empirical cumulative distribution function (CDF).

```R
R> plotdist(groundbeef$serving, histo = TRUE, demp = TRUE)
```

In addition to empirical plots, descriptive statistics may help to choose candidates to describe a distribution among a set of parametric distributions. Especially the skewness and kurtosis,
<figure><img src="image_5.png" /><figcaption>Figure 2: Skewness-kurtosis plot for a continuous variable (serving size from the groundbeef data set) as provided by the `descdist` function.</figcaption></figure>

linked to the third and fourth moments, are useful for this purpose. A non-zero skewness reveals a lack of symmetry of the empirical distribution, while the kurtosis value quantifies the weight of tails in comparison to the normal distribution for which the kurtosis equals 3. The skewness and kurtosis and their corresponding unbiased estimator (Casella and Berger 2002) from a sample $(X_i)_i \stackrel{\text{i.i.d.}}{\sim} X$ with observations $(x_i)_i$ are given by

$$ sk(X) = \frac{E[(X - E(X))^3]}{\mathrm{VAR}(X)^{\frac{3}{2}}}, \quad \widehat{sk} = \frac{\sqrt{n(n-1)}}{n-2} \times \frac{m_3}{m_2^{\frac{3}{2}}}, \qquad (1) $$

$$ kr(X) = \frac{E[(X - E(X))^4]}{\mathrm{VAR}(X)^2}, \quad \widehat{kr} = \frac{n-1}{(n-2)(n-3)} \left( (n+1) \times \frac{m_4}{m_2^2} - 3(n-1) \right) + 3, \qquad (2) $$

where $m_2, m_3, m_4$ denote empirical moments defined by $m_k = \frac{1}{n} \sum_{i=1}^{n} (x_i - \bar{x})^k$, with $x_i$ the $n$ observations of variable $x$ and $\bar{x}$ their mean value.

The `descdist` function provides classical descriptive statistics (minimum, maximum, median, mean, standard deviation), skewness and kurtosis. By default, unbiased estimations of the three last statistics are provided. Nevertheless, the argument `method` can be changed from `"unbiased"` (default) to `"sample"` to obtain them without correction for bias. A skewness-kurtosis plot such as the one proposed by Cullen and Frey (1999) is provided by the `descdist` function for the empirical distribution (see Figure 2 for the `groundbeef` data set). On this plot, values for common distributions are displayed in order to help the choice of distributions to fit to data. For some distributions (normal, uniform, logistic, exponential), there is only
one possible value for the skewness and the kurtosis. Thus, the distribution is represented by a single point on the plot. For other distributions, areas of possible values are represented, consisting in lines (as for the gamma and lognormal distributions), or larger areas (as for the beta distribution).

Skewness and kurtosis are known not to be robust. In order to take into account the uncertainty of the estimated values of kurtosis and skewness from data, a nonparametric bootstrap procedure (Efron and Tibshirani 1994) can be performed by using the argument boot. Values of skewness and kurtosis are computed on bootstrap samples (constructed by random sampling with replacement from the original data set) and reported on the skewness-kurtosis plot. Nevertheless, the user needs to know that skewness and kurtosis, like all higher moments, have a very high variance. This is a problem which cannot be completely solved by the use of bootstrap. The skewness-kurtosis plot should then be regarded as indicative only. The properties of the random variable should be considered, notably its expected value and its range, as a complement to the use of the plotdist and descdist functions. Below is a call to the descdist function to describe the distribution of the serving size from the groundbeef data set and to draw the corresponding skewness-kurtosis plot (see Figure 2). Looking at the results on this example with a positive skewness and a kurtosis not far from 3, the fit of three common right-skewed distributions could be considered, Weibull, gamma and lognormal distributions.

R> descdist(groundbeef$serving, boot = 1000)

summary statistics

------

min: 10 max: 200

median: 79

mean: 73.65

estimated sd: 35.88

estimated skewness: 0.7353

estimated kurtosis: 3.551

## 2.2. Fit of distributions by maximum likelihood estimation

Once selected, one or more parametric distributions $f(\cdot|\theta)$ (with parameter $\theta \in \mathbb{R}^d$) may be fitted to the data set, one at a time, using the fitdist function. Under the i.i.d. sample assumption, distribution parameters $\theta$ are by default estimated by maximizing the likelihood function defined as:

$$L(\theta) = \prod_{i=1}^{n} f(x_i|\theta) \qquad (3)$$

with $x_i$ the $n$ observations of variable $X$ and $f(\cdot|\theta)$ the density function of the parametric distribution. The other proposed estimation methods are described in Section 3.1.

The fitdist function returns an S3 object of class 'fitdist' for which print, summary and plot functions are provided. The fit of a distribution using fitdist assumes that the corresponding d, p, q functions (standing respectively for the density, the distribution and the quantile functions) are defined. Classical distributions are already defined in that way in the
stats package, e.g., dnorm, pnorm and qnorm for the normal distribution (see ?Distributions). Others may be found in various packages (see the CRAN Task View on Probability Distributions; Dutang 2014). Distributions not found in any package must be implemented by the user as d, p, q functions. In the call to fitdist, a distribution has to be specified via the argument dist either by the character string corresponding to its common root name used in the names of d, p, q functions (e.g., "norm" for the normal distribution) or by the density function itself, from which the root name is extracted (e.g., dnorm for the normal distribution). Numerical results returned by the fitdist function are (1) the parameter estimates, (2) the estimated standard errors (computed from the estimate of the Hessian matrix at the maximum likelihood solution), (3) the loglikelihood, (4) Akaike and Bayesian information criteria (the so-called AIC and BIC), and (5) the correlation matrix between parameter estimates. Below is a call to the fitdist function to fit a Weibull distribution to the serving size from the groundbeef data set.

R> fw <- fitdist(groundbeef$serving, "weibull")

R> summary(fw)

Fitting of the distribution ' weibull ' by maximum likelihood

Parameters :

estimate Std. Error

shape 2.186 0.1046

scale 83.348 2.5269

Loglikelihood: -1255 AIC: 2514 BIC: 2522

Correlation matrix:

shape scale

shape 1.0000 0.3218

scale 0.3218 1.0000

The plot of an object of class 'fitdist' provides four classical goodness-of-fit plots (Cullen and Frey 1999) presented in Figure 3:

• a density plot representing the density function of the fitted distribution along with the histogram of the empirical distribution,

• a CDF plot of both the empirical distribution and the fitted distribution,

• a Q-Q plot representing the empirical quantiles (y-axis) against the theoretical quantiles (x-axis)

• a P-P plot representing the empirical distribution function evaluated at each data point (y-axis) against the fitted distribution function (x-axis).

For CDF, Q-Q and P-P plots, the probability plotting position is defined by default using Hazen's rule, with probability points of the empirical distribution calculated as $(1:n - 0.5)/n$, as recommended by Blom (1959). This plotting position can be easily changed (see the reference manual for details; Delignette-Muller et al. 2015).

Unlike the generic plot function, the denscomp, cdfcomp, qqcomp and ppcomp functions enable to draw separately each of these four plots, in order to compare the empirical distribution and
<figure><img src="image_6.png" /><figcaption>Figure 3: Four goodness-of-fit plots for various distributions fitted to continuous data (Weibull, gamma and lognormal distributions fitted to serving sizes from the groundbeef data set) as provided by functions denscomp, qqcomp, cdfcomp and ppcomp.</figcaption></figure>

multiple parametric distributions fitted on the same data set. These functions must be called with the first argument corresponding to a list of objects of class 'fitdist', and optionally further arguments to customize the plot (see the reference manual for lists of arguments that may be specific to each plot; Delignette-Muller et al. 2015). In the following example, we compare the fit of a Weibull, a lognormal and a gamma distribution to the groundbeef data set (Figure 3).

```R
R> fg <- fitdist(groundbeef$serving, "gamma")
R> fln <- fitdist(groundbeef$serving, "lnorm")
R> par(mfrow = c(2, 2))
R> plot.legend <- c("Weibull", "lognormal", "gamma")
R> denscomp(list(fw, fln, fg), legendtext = plot.legend)
R> qqcomp(list(fw, fln, fg), legendtext = plot.legend)
R> cdfcomp(list(fw, fln, fg), legendtext = plot.legend)
R> ppcomp(list(fw, fln, fg), legendtext = plot.legend)
```
The density plot and the CDF plot may be considered as the basic classical goodness-of-fit plots. The two other plots are complementary and can be very informative in some cases. The Q-Q plot emphasizes the lack-of-fit at the distribution tails while the P-P plot emphasizes the lack-of-fit at the distribution center. In the present example (in Figure 3), none of the three fitted distributions correctly describes the center of the distribution, but the Weibull and gamma distributions could be preferred for their better description of the right tail of the empirical distribution, especially if this tail is important in the use of the fitted distribution, as it is in the context of food risk assessment.

The data set named endosulfan will now be used to illustrate other features of the fitdistrplus package. This data set contains acute toxicity values for the organochlorine pesticide endosulfan (geometric mean of LC50 or EC50 values in µg.L⁻<fnref n="1" />), tested on Australian and non-Australian laboratory-species (Hose and Van den Brink 2004). In ecotoxicology, a lognormal or a loglogistic distribution is often fitted to such a data set in order to characterize the species sensitivity distribution (SSD) for a pollutant. A low percentile of the fitted distribution, generally the 5% percentile, is then calculated and named the hazardous concentration 5% (HC5). It is interpreted as the value of the pollutant concentration protecting 95% of the species (Posthuma, Suter, and Traas 2010). But the fit of a lognormal or a loglogistic distribution to the whole endosulfan data set is rather bad (Figure 4), especially due to a minority of very high values. The two-parameter Pareto distribution and the three-parameter Burr distribution (which is an extension of both the loglogistic and the Pareto distributions) have been fitted. Pareto and Burr distributions are provided in the package actuar. So far, we did not have to define starting values (in the optimization process) as reasonable starting values are implicitly defined within the fitdist function for most of the distributions defined in R (see ?fitdist for details). For other distributions like the Pareto and the Burr distribution, initial values for the distribution parameters have to be supplied in the argument start, as a named list with initial values for each parameter (as they appear in the d, p, q functions). Having defined reasonable starting values<fnref n="1" />, various distributions can be fitted and graphically compared. In this example, the function cdfcomp can be used to report CDF values in a logscale so as to emphasize discrepancies on the tail of interest while defining an HC5 value (Figure 4).

```R
R> data("endosulfan", package = "fitdistrplus")
R> ATV <- endosulfan$ATV
R> fendo.ln <- fitdist(ATV, "lnorm")
R> library("actuar")
R> fendo.ll <- fitdist(ATV, "llogis", start = list(shape = 1, scale = 500))
R> fendo.P <- fitdist(ATV, "pareto", start = list.shape = 1, scale = 500))
R> fendo.B <- fitdist(ATV, "burr", start = list.shape1 = 0.3, shape2 = 1,
+ rate = 1))
R> cdfcomp(list(fendo.ln, fendo.ll, fendo.P, fendo.B), xlogscale = TRUE,
+ ylogscale = TRUE, legendtext = c("lognormal", "loglogistic", "Pareto",
+ "Burr"))
```

None of the fitted distributions correctly describes the right tail observed in the data set, but

<footnote n="1">¹The plotdist function can plot any parametric distribution with specified parameter values in argument para. It can thus help to find correct initial values for the distribution parameters in non trivial cases, by iterative calls if necessary (see the reference manual for examples; Delignette-Muller et al. 2015).</footnote>
<figure><img src="image_7.png" /><figcaption>Figure 4: CDF plot to compare the fit of four distributions to acute toxicity values of various organisms for the organochlorine pesticide endosulfan (endosulfan data set) as provided by the cdfcomp function, with CDF values in a logscale to emphasize discrepancies on the left tail.</figcaption></figure>

as shown in Figure 4, the left-tail seems to be better described by the Burr distribution. Its use could then be considered to estimate the HC5 value as the 5% quantile of the distribution. This can be easily done using the `quantile` generic function defined for an object of class 'fitdist'. Below is this calculation together with the calculation of the empirical quantile for comparison.

R> quantile(fendo.B, prob = 0.05)

Estimated quantiles for each specified probability (non-censored data)

p=0.05

estimate 0.2939

R> quantile(ATV, prob = 0.05)

5%

0.2

In addition to the ecotoxicology context, the `quantile` generic function is also attractive in the actuarial-financial context. In fact, the value-at-risk $VAR_α$ is defined as the $(1 - α)$-quantile of the loss distribution and can be computed with `quantile` on a 'fitdist' object.

The computation of different goodness-of-fit statistics is proposed in the `fitdistrplus` package in order to further compare fitted distributions. The purpose of goodness-of-fit statistics aims at measuring the distance between the fitted parametric distribution and the empirical distribution: e.g., the distance between the fitted cumulative distribution function $F$ and the empirical distribution function $F_n$. When fitting continuous distributions, three
<table><thead><tr><th>Statistic</th><th>General formula</th><th>Computational formula</th></tr></thead><tbody><tr><td>Kolmogorov-Smirnov<br>(KS)</td><td>sup |F<sub>n</sub>(x) - F(x)|</td><td>max(D<sup>+</sup>, D<sup>-</sup>) with<br>D<sup>+</sup> = max<sub>i=1,...,n</sub> (<sup>i</sup>/<sub>n</sub> - F<sub>i</sub>)<br>D<sup>-</sup> = max<sub>i=1,...,n</sub> (F<sub>i</sub> - <sup>i-1</sup>/<sub>n</sub>)</td></tr><tr><td>Cramer-von Mises<br>(CvM)</td><td>n &int;<sub>-&infin;</sub><sup>&infin;</sup> (F<sub>n</sub>(x) - F(x))<sup>2</sup> dx</td><td><sup>1</sup>/<sub>12n</sub> + &sum;<sub>i=1</sub><sup>n</sup> (F<sub>i</sub> - <sup>2i-1</sup>/<sub>2n</sub>)<sup>2</sup></td></tr><tr><td>Anderson-Darling<br>(AD)</td><td>n &int;<sub>-&infin;</sub><sup>&infin;</sup> (F<sub>n</sub>(x)-F(x))<sup>2</sup> / (F(x)(1-F(x))) dx</td><td>-n - <sup>1</sup>/<sub>n</sub> &sum;<sub>i=1</sub><sup>n</sup> (2i - 1) log(F<sub>i</sub>(1-F<sub>n+1-i</sub>))</td></tr></tbody></table>

where $F_i \triangleq F(x_i)$

<caption>Table 1: Goodness-of-fit statistics as defined by D'Agostino and Stephens (1986).</caption>

goodness-of-fit statistics are classically considered: Cramer-von Mises, Kolmogorov-Smirnov and Anderson-Darling statistics (D'Agostino and Stephens 1986). Naming $x_i$ the n observations of a continuous variable X arranged in an ascending order, Table 1 gives the definition and the empirical estimate of the three considered goodness-of-fit statistics. They can be computed using the function `gofstat` as defined by D'Agostino and Stephens (1986).

R> `gofstat(list(fendo.ln, fendo.ll, fendo.P, fendo.B),` + `fitnames = c("lnnorm", "llogis", "Pareto", "Burr"))`

Goodness-of-fit statistics

<table><thead><tr><th></th><th>lnnorm</th><th>llogis</th><th>Pareto</th><th>Burr</th></tr></thead><tbody><tr><td>Kolmogorov-Smirnov statistic</td><td>0.1672</td><td>0.1196</td><td>0.08488</td><td>0.06155</td></tr><tr><td>Cramer-von Mises statistic</td><td>0.6374</td><td>0.3827</td><td>0.13926</td><td>0.06803</td></tr><tr><td>Anderson-Darling statistic</td><td>3.4721</td><td>2.8316</td><td>0.89206</td><td>0.52393</td></tr></tbody></table>

Goodness-of-fit criteria

<table><thead><tr><th></th><th>lnnorm</th><th>llogis</th><th>Pareto</th><th>Burr</th></tr></thead><tbody><tr><td>Aikake's Information Criterion</td><td>1069</td><td>1069</td><td>1048</td><td>1046</td></tr><tr><td>Bayesian Information Criterion</td><td>1074</td><td>1075</td><td>1053</td><td>1054</td></tr></tbody></table>

Due to giving more weight to distribution tails, the Anderson-Darling statistic is of special interest when it matters to equally emphasize the tails as well as the main body of a distribution. This is often the case in risk assessment (Cullen and Frey 1999; Vose 2010). For this reason, this statistic is often used to select the best distribution among those fitted. Nevertheless, this statistic should be used cautiously when comparing fits of various distributions. Keeping in mind that the weighting of each CDF quadratic difference depends on the parametric distribution in its definition (see Table 1), Anderson-Darling statistics computed for several distributions fitted on the same data set are theoretically difficult to compare. Moreover, such a statistic, as the Cramer-von Mises and Kolmogorov-Smirnov ones, does not take into account the complexity of the model (i.e., parameter number). It is not a problem when the compared distributions are characterized by the same number of parameters, but it could systematically promote the selection of the more complex distributions otherwise. Looking at
classical penalized criteria based on the loglikelihood (AIC, BIC) seems thus also interesting, especially to discourage overfitting.

In the previous example, all the goodness-of-fit statistics based on the CDF distance are in favor of the Burr distribution, the only one characterized by three parameters, while AIC and BIC values respectively give the preference to the Burr distribution or the Pareto distribution. The choice between these two distributions seems thus less obvious and could be discussed. Even if specifically recommended for discrete distributions, the Chi-squared statistic may also be used for continuous distributions (see Section 3.3 and the reference manual for examples; Delignette-Muller *et al.* 2015).

## 2.3. Uncertainty in parameter estimates

The uncertainty in the parameters of the fitted distribution can be estimated by parametric or nonparametric bootstraps using the `boodist` function for non-censored data (Efron and Tibshirani 1994). This function returns the bootstrapped values of parameters in an S3 class object which can be plotted to visualize the bootstrap region. The medians and the 95 percent confidence intervals of parameters (2.5 and 97.5 percentiles) are printed in the summary. When smaller than the total number of iterations (due to lack of convergence of the optimization algorithm for some bootstrapped data sets), the number of iterations for which the estimation converges is also printed in the summary.

The plot of an object of class 'bootdist' consists in a scatter plot or a matrix of scatter plots of the bootstrapped values of parameters providing a representation of the joint uncertainty distribution of the fitted parameters. Below is an example of the use of the bootdist function with the previous fit of the Burr distribution to the endosulfan data set (Figure 5).

R> bendo.B <- bootdist(fendo.B, niter = 1001)

R> summary(bendo.B)

Parametric bootstrap medians and 95% percentile CI

Median 2.5% 97.5%

shape1 0.2008 0.09968 0.3625

shape2 1.5803 1.06174 2.7959

rate 1.4859 0.69571 2.7430

The estimation method converged only for 1000 among 1001 iterations

R> plot(bendo.B)

Bootstrap samples of parameter estimates are useful especially to calculate confidence intervals on each parameter of the fitted distribution from the marginal distribution of the bootstrapped values. It is also interesting to look at the joint distribution of the bootstrapped values in a scatter plot (or a matrix of scatter plots if the number of parameters exceeds two) in order to understand the potential structural correlation between parameters (see Figure 5).

The use of the whole bootstrap sample is also of interest in the risk assessment field. Its use enables the characterization of uncertainty in distribution parameters. It can be directly used within a second-order Monte Carlo simulation framework, especially within the package `mc2d`
<figure><img src="image_8.png" /><figcaption>Figure 5: Bootstrapped values of parameters for a fit of the Burr distribution characterized by three parameters (example on the *endosulfan* data set) as provided by the plot of an object of class 'bootdist'.</figcaption></figure>

(Pouillot, Delignette-Muller, and Denis 2011). One could refer to Pouillot and Delignette-Muller (2010) for an introduction to the use of **mc2d** and **fitdistrplus** packages in the context of quantitative risk assessment.

The bootstrap method can also be used to calculate confidence intervals on quantiles of the fitted distribution. For this purpose, a generic **quantile** function is provided for class 'bootdist'. By default, 95% percentile bootstrap confidence intervals of quantiles are provided. Going back to the previous example from ecotoxicology, this function can be used to estimate the uncertainty associated to the HC5 estimation, for example from the previously fitted Burr distribution to the *endosulfan* data set.

R> quantile(bendo.B, prob = 0.05)

(original) estimated quantiles for each specified probability (non-censored data)

p=0.05

estimate 0.2939

Median of bootstrap estimates

p=0.05

estimate 0.3004

two-sided 95 % CI of each quantile
p=0.05
2.5 % 0.1785
97.5 % 0.5000

The estimation method converged only for 1000 among 1001 bootstrap iterations.

### 3. Advanced topics

#### 3.1. Alternative methods for parameter estimation

This subsection focuses on alternative estimation methods. One of the alternatives for continuous distributions is the maximum goodness-of-fit estimation method also called minimum distance estimation method (D'Agostino and Stephens 1986; Dutang et al. 2008). In this package this method is proposed with eight different distances: the three classical distances defined in Table 1, or one of the variants of the Anderson-Darling distance proposed by Luceno (2006) and defined in Table 2. The right-tail Anderson-Darling (AD) gives more weight to the right-tail, the left-tail AD gives more weight to the left tail. Either of the tails, or both of them, can receive even larger weights by using second order AD statistics.

To fit a distribution by maximum goodness-of-fit estimation, one needs to fix the argument method to "mge" in the call to fitdist and to specify the argument gof coding for the chosen goodness-of-fit distance. This function is intended to be used only with continuous non-censored data.

Maximum goodness-of-fit estimation may be useful to give more weight to data at one tail of the distribution. In the previous example from ecotoxicology, we used a non classical

<table><thead><tr><th>Statistic</th><th>General formula</th><th>Computational formula</th></tr></thead><tbody><tr><td>Right-tail AD<br>(ADR)</td><td>∫<sub>-∞</sub><sup>∞</sup> <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 0.2em;">(F<sub>n</sub>(x)-F(x))<sup>2</sup></div><div style="padding: 0 0.2em;">1-F(x)</div></div>dx</td><td><div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 0.2em;">n</div><div style="padding: 0 0.2em;">2</div></div> - 2 ∑<sub>i=1</sub><sup>n</sup> F<sub>i</sub> - <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 0.2em;">1</div><div style="padding: 0 0.2em;">n</div></div> ∑<sub>i=1</sub><sup>n</sup> (2i - 1) ln(F&#773;<sub>n+1-i</sub>)</td></tr><tr><td>Left-tail AD<br>(ADL)</td><td>∫<sub>-∞</sub><sup>∞</sup> <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 0.2em;">(F<sub>n</sub>(x)-F(x))<sup>2</sup></div><div style="padding: 0 0.2em;">(F(x))</div></div>dx</td><td>- <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 0.2em;">3n</div><div style="padding: 0 0.2em;">2</div></div> + 2 ∑<sub>i=1</sub><sup>n</sup> F<sub>i</sub> - <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 0.2em;">1</div><div style="padding: 0 0.2em;">n</div></div> ∑<sub>i=1</sub><sup>n</sup> (2i - 1) ln(F<sub>i</sub>)</td></tr><tr><td>Right-tail AD<br>2nd order (AD2R)</td><td>ad2r = ∫<sub>-∞</sub><sup>∞</sup> <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 0.2em;">(F<sub>n</sub>(x)-F(x))<sup>2</sup></div><div style="padding: 0 0.2em;">(1-F(x))<sup>2</sup></div></div>dx</td><td>ad2r = 2 ∑<sub>i=1</sub><sup>n</sup> ln(F&#773;<sub>i</sub>) + <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 0.2em;">1</div><div style="padding: 0 0.2em;">n</div></div> ∑<sub>i=1</sub><sup>n</sup> <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 0.2em;">2i-1</div><div style="padding: 0 0.2em;">F&#773;<sub>n+1-i</sub></div></div></td></tr><tr><td>Left-tail AD<br>2nd order (AD2L)</td><td>ad2l = ∫<sub>-∞</sub><sup>∞</sup> <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 0.2em;">(F<sub>n</sub>(x)-F(x))<sup>2</sup></div><div style="padding: 0 0.2em;">(F(x))<sup>2</sup></div></div>dx</td><td>ad2l = 2 ∑<sub>i=1</sub><sup>n</sup> ln(F<sub>i</sub>) + <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 0.2em;">1</div><div style="padding: 0 0.2em;">n</div></div> ∑<sub>i=1</sub><sup>n</sup> <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 0.2em;">2i-1</div><div style="padding: 0 0.2em;">F<sub>i</sub></div></div></td></tr><tr><td>AD 2nd order<br>(AD2)</td><td>ad2r + ad2l</td><td>ad2r + ad2l</td></tr></tbody></table>

where $F_i \triangleq F(x_i)$; $\bar{F}_i \triangleq 1 - F(x_i)$

<caption>Table 2: Modified Anderson-Darling statistics as defined by Luceno (2006).</caption>
<figure><img src="image_9.png" /><figcaption>Figure 6: Comparison of a lognormal distribution fitted by MLE and by MGE using two different goodness-of-fit distances: left-tail Anderson-Darling and left-tail Anderson Darling of second order (example with the endosulfan data set) as provided by the cdfcomp function, with CDF values in a logscale to emphasize discrepancies on the left tail.</figcaption></figure>

distribution (the Burr distribution) to correctly fit the empirical distribution especially on its left tail. In order to correctly estimate the 5% percentile, we could also consider the fit of the classical lognormal distribution, but minimizing a goodness-of-fit distance giving more weight to the left tail of the empirical distribution. In what follows, the left tail AD distances of first or second order are used to fit a lognormal to endosulfan data set (see Figure 6).

```R
R> fendo.ln.ADL <- fitdist(ATV, "lnnorm", method = "mge", gof = "ADL")
R> fendo.ln.AD2L <- fitdist(ATV, "lnnorm", method = "mge", gof = "AD2L")
R> cdfcomp(list(fendo.ln, fendo.ln.ADL, fendo.ln.AD2L),
+           xlogscale = TRUE, ylogscale = TRUE,
+           main = "Fitting a lognormal distribution", xlegend = "bottomright",
+           legendtext = c("MLE", "Left-tail AD", "Left-tail AD 2nd order"))
```

Comparing the 5% percentiles (HC5) calculated using these three fits to the one calculated from the MLE fit of the Burr distribution, we can observe, on this example, that fitting the lognormal distribution by maximizing left tail Anderson-Darling distances of first or second order enables to approach the value obtained by fitting the Burr distribution by MLE.

```R
R> (HC5.estimates <- c( empirical = as.numeric(quantile(ATV, EV, EV, 0.05)),
+                       Burr = as.numeric(quantile(fendo.B, EV, EV, 0.05)$quantiles),
+                       lognormal_MLE = as.numeric(quantile(fendo.ln, EV, EV, 0.05)$quantiles),
+                       lognormal_AD2 = as.numeric(quantile(fendo.ln.ADL,
+                                         EV, EV, 0.05)$quantiles),
+                       lognormal_AD2L = as.numeric(quantile(fendo.ln.AD2L,
+                                         EV, EV, 0.05)$quantiles)))
```
<table><thead><tr><th>empirical</th><th>Burr</th><th>lognormal_MLE</th><th>lognormal_AD2</th><th>lognormal_AD2L</th></tr></thead><tbody><tr><td>0.20000</td><td>0.29393</td><td>0.07259</td><td>0.19591</td><td>0.25877</td></tr></tbody></table>

The moment matching estimation (MME) is another method commonly used to fit parametric distributions (Vose 2010). MME consists in finding the value of the parameter $\theta$ that equalizes the first theoretical raw moments of the parametric distribution to the corresponding empirical raw moments as in Equation 4:

$$E(X^k|\theta) = \frac{1}{n} \sum_{i=1}^{n} x_i^k, \qquad (4)$$

for $k = 1, \dots, d$, with $d$ the number of parameters to estimate and $x_i$ the $n$ observations of variable $X$. For moments of order greater than or equal to 2, it may also be relevant to match centered moments. Therefore, we match the moments given in Equation 5:

$$E(X|\theta) = \bar{x}, \quad E((X - E(X))^k|\theta) = m_k, \quad \text{for } k = 2, \dots, d, \qquad (5)$$

where $m_k$ denotes the empirical centered moments. This method can be performed by setting the argument `method` to `"mme"` in the call to `fitdist`. The estimate is computed by a closed-form formula for the following distributions: normal, lognormal, exponential, Poisson, gamma, logistic, negative binomial, geometric, beta and uniform distributions. In this case, for distributions characterized by one parameter (geometric, Poisson and exponential), this parameter is simply estimated by matching theoretical and observed means, and for distributions characterized by two parameters, these parameters are estimated by matching theoretical and observed means and variances (Vose 2010). For other distributions, the equation of moments is solved numerically using the `optim` function by minimizing the sum of squared differences between observed and theoretical moments (see the `fitdistrplus` reference manual for technical details; Delignette-Muller et al. 2015).

A classical data set from the Danish insurance industry published in McNeil (1997) will be used to illustrate this method. In `fitdistrplus`, the data set is stored in `danishuni` for the univariate version and contains the loss amounts collected at Copenhagen Reinsurance between 1980 and 1990. In actuarial science, it is standard to consider positive heavy-tailed distributions and have a special focus on the right-tail of the distributions. In this numerical experiment, we choose classic actuarial distributions for loss modeling: the lognormal distribution and the Pareto type II distribution (Klugman, Panjer, and Willmot 2009).

The lognormal distribution is fitted to `danishuni` data set by matching moments implemented as a closed-form formula. On the left-hand graph of Figure 7, the fitted distribution functions obtained using the moment matching estimation (MME) and maximum likelihood estimation (MLE) methods are compared. The MME method provides a more cautious estimation of the insurance risk as the MME-fitted distribution function (resp. MLE-fitted) underestimates (overestimates) the empirical distribution function for large values of claim amounts.

R> data("danishuni", package = "fitdistrplus")

R> str(danishuni)

'data.frame': 2167 obs. of 2 variables:
$ Date: Date, format: "1980-01-03" "1980-01-04" ...
$ Loss: num 1.68 2.09 1.73 1.78 4.61 ...
<figure><img src="image_10.png" /><figcaption>Figure 7: Comparison between MME and MLE when fitting a lognormal or a Pareto distribution to loss data from the danishuni data set.</figcaption></figure>

```R
R> fdanish.ln.MLE <- fitdist(danishuni$Loss, "lnnorm")
R> fdanish.ln.MME <- fitdist(danishuni$Loss, "lnnorm", method = "mme",
+     order = 1:2)
R> cdfcomp(list(fdanish.ln.MLE, fdanish.ln.MME), legend = c("lognormal MLE",
+ "lognormal MME"), main = "Fitting a lognormal distribution",
+ xlogscale = TRUE, datapch = 20)
```

Next, a Pareto distribution, which gives more weight to the right-tail of the distribution, is fitted. As the lognormal distribution, the Pareto has two parameters, which allows a fair comparison.

We use the implementation of the `actuar` package providing raw and centered moments for that distribution (in addition to `d`, `p`, `q` and `r` functions; Goulet 2014). Fitting a heavy-tailed distribution for which the first and the second moments do not exist for certain values of the shape parameter requires some cautiousness. This is carried out by providing, for the optimization process, a lower and an upper bound for each parameter. The code below calls the L-BFGS-B optimization method in `optim`, since this quasi-Newton method allows for box constraints<fnref n="2" />. We choose to match moments defined by Equation 4, and so a function for computing the empirical raw moment (called `memp` in our example) is passed to `fitdist`. For two-parameter distributions (i.e., $d = 2$), Equations 4 and 5 are equivalent.

```R
R> library("actuar")
R> fdanish.P.MLE <- fitdist(danishuni$Loss, "pareto",
+     start = c(shape = 10, scale = 10), lower = 2+1e-6, upper = Inf)
R> memp <- function(x, order) sum(x^order)/length(x)
R> fdanish.P.MME <- fitdist(danishuni$Loss, "pareto", method = "mme",
+     order = 1:2, memp = "memp", start = c.shape = 10, scale = 10),
+     lower = c(2+1e-6, 2+1e-6), upper = c(Inf, Inf))
```

<footnote n="2">²That is what the B stands for.</footnote>
```R
R> cdfcomp(list(fdanish.P.MLE, fdanish.P.MME), legend = c("Pareto MLE",
+ "Pareto MME"), main = "Fitting a Pareto distribution",
+ xlogscale = TRUE, datapch = ".")
R> gofstat(list(fdanish.ln.MLE, fdanish.P.MLE, fdanish.ln.MME,
+ fdanish.P.MME), fitnames = c("lnnorm.mle", "Pareto.mle", "lnnorm.mme",
+ "Pareto.mme"))
```

Goodness-of-fit statistics

<table><thead><tr><th></th><th>lnnorm.mle</th><th>Pareto.mle</th><th>lnnorm.mme</th><th>Pareto.mme</th></tr></thead><tbody><tr><td>Kolmogorov-Smirnov statistic</td><td>0.1375</td><td>0.3124</td><td>0.4368</td><td>0.3638</td></tr><tr><td>Cramer-von Mises statistic</td><td>14.7911</td><td>37.7166</td><td>88.9503</td><td>53.0783</td></tr><tr><td>Anderson-Darling statistic</td><td>87.1933</td><td>208.3139</td><td>416.2567</td><td>272.4729</td></tr></tbody></table>

Goodness-of-fit criteria

<table><thead><tr><th></th><th>lnnorm.mle</th><th>Pareto.mle</th><th>lnnorm.mme</th><th>Pareto.mme</th></tr></thead><tbody><tr><td>Aikake's Information Criterion</td><td>8120</td><td>9250</td><td>9792</td><td>9395</td></tr><tr><td>Bayesian Information Criterion</td><td>8131</td><td>9261</td><td>9803</td><td>9407</td></tr></tbody></table>

As shown in Figure 7, MME and MLE fits are far less distant (when looking at the right-tail) for the Pareto distribution than for the lognormal distribution on this data set. Furthermore, for these two distributions, the MME method better fits the right-tail of the distribution from a visual point of view. This seems logical since empirical moments are influenced by large observed values. In the previous traces, we gave the values of goodness-of-fit statistics. Whatever the statistic considered, the MLE-fitted lognormal always provides the best fit to the observed data.

Maximum likelihood and moment matching estimations are certainly the most commonly used methods for fitting distributions (Cullen and Frey 1999). Keeping in mind that these two methods may produce very different results, the user should be aware of its great sensitivity to outliers when choosing the moment matching estimation. This may be seen as an advantage in our example if the objective is to better describe the right tail of the distribution, but it may be seen as a drawback if the objective is different.

Fitting of a parametric distribution may also be done by matching theoretical quantiles of the parametric distributions (for specified probabilities) against the empirical quantiles (Tse 2009). The equality of theoretical and empirical quantiles is expressed by Equation 6 below, which is very similar to Equations 4 and 5:

$$F^{-1}(p_k|\theta) = Q_{n,p_k} \qquad (6)$$

for $k = 1, \dots, d$, with $d$ the number of parameters to estimate (dimension of $\theta$ if there is no fixed parameters) and $Q_{n,p_k}$ the empirical quantiles calculated from data for specified probabilities $p_k$.

Quantile matching estimation (QME) is performed by setting the argument `method` to `"qme"` in the call to `fitdist` and adding an argument `probs` defining the probabilities for which the quantile matching is performed. The length of this vector must be equal to the number of parameters to estimate (as the vector of moment orders for MME). Empirical quantiles are computed using the `quantile` function of the `stats` package using `type = 7` by default (see `?quantile` and Hyndman and Fan 1996). But the type of quantile can be easily changed by
<figure><img src="image_11.png" /><figcaption>Figure 8: Comparison between QME and MLE when fitting a lognormal distribution to loss data from the danishuni data set.</figcaption></figure>

using the `qty` argument in the call to the `qme` function. The quantile matching is carried out numerically, by minimizing the sum of squared differences between observed and theoretical quantiles.

```r
R> fdanish.ln.QME1 <- fitdist(danishuni$Loss, "lnorm", method = "qme",
+   prob = c(1/3, 2/3))
R> fdanish.ln.QME2 <- fitdist(danishuni$Loss, "lnorm", method = "qme",
+   prob = c(8/10, 9/10))
R> cdfcomp(list(fdanish.ln.MLE, fdanish.ln.QME1, fdanish.ln.QME2),
+   legend = c("MLE", "QME(1/3, 2/3)", "QME(8/10, 9/10)"),
+   main = "Fitting a lognormal distribution", xlogscale = TRUE,
+   datapch = 20)
```

Above is an example of fitting a lognormal distribution to the danishuni data set by matching probabilities ($p_1 = 1/3, p_2 = 2/3$) and ($p_1 = 8/10, p_2 = 9/10$). As expected, the second QME fit gives more weight to the right-tail of the distribution. Compared to the maximum likelihood estimation, the second QME fit best suits the right-tail of the distribution, whereas the first QME fit best models the body of the distribution. The quantile matching estimation is of particular interest when we need to focus around particular quantiles, e.g., $p = 99.5\%$ in the Solvency II insurance context or $p = 5\%$ for the HC5 estimation in the ecotoxicology context.

## 3.2. Customization of the optimization algorithm

Each time a numerical minimization is carried out in the `fitdistrplus` package, the `optim` function of the `stats` package is used by default with the "Nelder-Mead" method for distributions characterized by more than one parameter and the "BFGS" method for distributions characterized by only one parameter. Sometimes the default algorithm fails to converge. It is then interesting to change some options of the `optim` function or to use another optimization func-
tion than optim to minimize the objective function. The argument optim.method can be used in the call to fitdist or fitdistcens. It will internally be passed to mledist, mmedist, mgedist or qmedist, and to optim (see ?optim for details about the different algorithms available).

Even if no error is raised when computing the optimization, changing the algorithm is of particular interest to enforce bounds on some parameters. For instance, a volatility parameter $\sigma$ is strictly positive $\sigma > 0$ and a probability parameter $p$ lies in $p \in [0, 1]$. This is possible by using arguments lower and/or upper, which – if used – automatically force optim.method = "L-BFGS-B".

Below are examples of fits of a gamma distribution $\mathcal{G}(\alpha, \lambda)$ to the groundbeef data set with various algorithms. Note that the conjugate gradient algorithm ("CG") needs far more iterations to converge (around 2500 iterations) compared to other algorithms (converging in less than 100 iterations).

```R
R> data("groundbeef", package = "fitdistrplus")
R> fNM <- fitdist(groundbeef$serving, "gamma", optim.method = "Nelder-Mead")
R> fBFGS <- fitdist(groundbeef$serving, "gamma", optim.method = "BFGS")
R> fSANN <- fitdist(groundbeef$serving, "gamma", optim.method = "SANN")
R> fCG <- try(fitdist(groundbeef$serving, "gamma", optim.method = "CG"),
+           control = list(maxit = 10000)))
R> if (is(fCG, "try-error"))
+     fCG <- list(estimate = NA)
```

It is also possible to use another function than optim to minimize the objective function by specifying the argument custom.optim in the call to fitdist. It may be necessary to customize this optimization function to meet the following requirements. (1) the custom.optim function must have the following arguments: fn for the function to be optimized and par for the initialized parameters. (2) custom.optim should carry out a minimization and must return the following components: par for the estimate, convergence for the convergence code, value = fn(par) and hessian. Below is an example of code written to wrap the genoud function from the rgenoud package (Mebane, Jr. and Sekhon 2011) in order to respect our optimization “template”. The rgenoud package implements the genetic (stochastic) algorithm.

```R
R> mygenoud <- function(fn, par, ...) {
+     require("rgenoud")
+     res <- genoud(fn, starting.values = par, ...)
+     standardres <- c(res, convergence = 0)
+     return(standardres)
+ }
```

The customized optimization function can then be passed as the argument custom.optim in the call to fitdist or fitdistcens. The following code can for example be used to fit a gamma distribution to the groundbeef data set. Note that in this example various arguments are also passed from fitdist to genoud: nvars, Domains, boundary.enforcement, print.level and hessian. The code below compares all the parameter estimates ($\hat{\alpha}, \hat{\lambda}$) by the different algorithms: shape $\alpha$ and rate $\lambda$ parameters are relatively similar on this example, roughly 4.00 and 0.05, respectively.
```R
R> fgenoud <- mledist(groundbeef$serving, "gamma", custom.optim = mygenoud,
+                      nvars = 2, max.generations = 10, Domains = cbind(c(0, 0), c(10, 10)),
+                      boundary.enforcement = 1, hessian = TRUE, print.level = 0, P9 = 10)
R> cbind(NM = fNM$estimate, BFGS = fBFGS$estimate, SANN = fSANN$estimate,
+       CG = fCG$estimate, fgenoud = fgenoud$estimate)
```

```
NM      BFGS     SANN      CG fgenoud
shape  4.00825  4.22848  3.97131  4.12850  4.00834
rate   0.05442  0.05742  0.05392  0.05606  0.05443
```

### 3.3. Fitting distributions to other types of data

Analytical methods often lead to semi-quantitative results which are referred to as censored data. Observations only known to be under a limit of detection are left-censored data. Observations only known to be above a limit of quantification are right-censored data. Results known to lie between two bounds are interval-censored data. These two bounds may correspond to a limit of detection and a limit of quantification, or more generally to uncertainty bounds around the observation. Right-censored data are also commonly encountered with survival data (Klein and Moeschberger 2003). A data set may thus contain right-, left-, or interval-censored data, or may be a mixture of these categories, possibly with different upper and lower bounds. Censored data are sometimes excluded from the data analysis or replaced by a fixed value, which in both cases may lead to biased results. A more recommended approach to correctly model such data is based upon maximum likelihood (Klein and Moeschberger 2003; Helsel 2005).

Censored data may thus contain left-censored, right-censored and interval-censored values, with several lower and upper bounds. Before their use in package **fitdistrplus**, such data must be coded into a dataframe with two columns, respectively named **left** and **right**, describing each observed value as an interval. The **left** column contains either **NA** for left censored observations, the left bound of the interval for interval censored observations, or the observed value for non-censored observations. The **right** column contains either **NA** for right censored observations, the right bound of the interval for interval censored observations, or the observed value for non-censored observations. To illustrate the use of package **fitdistrplus** to fit distributions to censored continuous data, we will use another data set from ecotoxicology, included in our package and named **salinity**. This data set contains acute salinity tolerance (LC50 values in electrical conductivity, $mS.cm^{-1}$) of riverine macro-invertebrates taxa from the southern Murray-Darling Basin in Central Victoria, Australia (Kefford, Fields, Clay, and Nugegoda 2007).

```R
R> data("salinity", package = "fitdistrplus")
R> str(salinity)
```

```
'data.frame':      108 obs. of 2 variables:
$ left : num  20 20 20 20 20 21.5 15 20 23.7 25 ...
$ right: num  NA NA NA NA NA 21.5 30 25 23.7 NA ...
```

Using censored data such as those coded in the **salinity** data set, the empirical distribution can be plotted using the **plotdistcens** function. By default, this function uses
<figure><img src="image_12.png" /><figcaption>Figure 9: Simple plot of censored raw data (72-hour acute salinity tolerance of riverine macro-invertebrates from the salinity data set) as ordered points and intervals.</figcaption></figure>

the expectation-maximization approach of Turnbull (1974) to compute the overall empiri-cal CDF curve with optional confidence intervals, by calls to `survfit` and the `plot` method for 'survfit' objects from the `survival` package (Therneau 2014; Therneau and Grambsch 2000, Figure 10 shows the Turnbull plot of data together with two fitted distributions). A less rigorous but sometimes more illustrative plot can be obtained by fixing the argument `Turnbull` to `FALSE` in the call to `plotdistcens` (see Figure 9 for an example and the help page of Function `plotdistcens` for details). This plot enables to see the real nature of censored data, as points and intervals.

R> plotdistcens(salinity, Turnbull = FALSE)

As for non censored data, one or more parametric distributions can be fitted to the censored data set, one at a time, but using in this case the `fitdistcens` function. This function estimates the vector of distribution parameters $\theta$ by maximizing the likelihood for censored data defined as:

$$L(\theta) = \prod_{i=1}^{N_{nonC}} f(x_i|\theta) \times \prod_{j=1}^{N_{leftC}} F(x_j^{upper}|\theta) \times \prod_{k=1}^{N_{rightC}} (1 - F(x_k^{lower}|\theta)) \times \prod_{m=1}^{N_{intC}} (F(x_m^{upper}|\theta) - F(x_m^{lower}|\theta)) \quad (7)$$

with $x_i$ the $N_{nonC}$ non-censored observations, $x_j^{upper}$ upper values defining the $N_{leftC}$ left-censored observations, $x_k^{lower}$ lower values defining the $N_{rightC}$ right-censored observations, $[x_m^{lower}, x_m^{upper}]$ the intervals defining the $N_{intC}$ interval-censored observations, and $F$ the cumulative distribution function of the parametric distribution (Klein and Moeschberger 2003; Helsel 2005).
As `fitdist`, `fitdistcens` returns the results of the fit of any parametric distribution to a data set as an S3 class object that can be easily printed, summarized or plotted. For the `salinity` data set, a lognormal distribution or a loglogistic can be fitted as commonly done in ecotoxicology for such data. As with `fitdist`, for some distributions (see Delignette-Muller et al. 2015, for details), it is necessary to specify initial values for the distribution parameters in the argument `start`. The `plotdistcens` function can help to find correct initial values for the distribution parameters in non trivial cases, by a manual iterative use if necessary.

R> fsal.ln <- fitdistcens(salinity, "lnnorm")

R> fsal.l1 <- fitdistcens(salinity, "llogis", start = list shape = 5,
+ scale = 40))

R> summary(fsal.ln)

FITTING OF THE DISTRIBUTION ' lnorm ' BY MAXIMUM LIKELIHOOD ON CENSORED DATA
PARAMETERS

<table><thead><tr><th></th><th>estimate</th><th>Std. Error</th></tr></thead><tbody><tr><td>meanlog</td><td>3.3854</td><td>0.06487</td></tr><tr><td>sdlog</td><td>0.4961</td><td>0.05455</td></tr></tbody><tfoot><tr><td colspan="3">Loglikelihood: -139.1 AIC: 282.1 BIC: 287.5</td></tr></tfoot></table>

Correlation matrix:

<table><tr><td>meanlog</td><td>sdlog</td></tr><tr><td>meanlog</td><td>1.0000</td><td>0.2938</td></tr><tr><td>sdlog</td><td>0.2938</td><td>1.0000</td></tr></table>

R> summary(fsal.l1)

FITTING OF THE DISTRIBUTION ' llogis ' BY MAXIMUM LIKELIHOOD ON CENSORED DATA
PARAMETERS

<table><thead><tr><th></th><th>estimate</th><th>Std. Error</th></tr></thead><tbody><tr><td>shape</td><td>3.421</td><td>0.4158</td></tr><tr><td>scale</td><td>29.930</td><td>1.9447</td></tr></tbody><tfoot><tr><td colspan="3">Loglikelihood: -140.1 AIC: 284.1 BIC: 289.5</td></tr></tfoot></table>

Correlation matrix:

<table><tr><td>shape</td><td>scale</td></tr><tr><td>shape</td><td>1.0000</td><td>-0.2022</td></tr><tr><td>scale</td><td>-0.2022</td><td>1.0000</td></tr></table>

Computations of goodness-of-fit statistics have not yet been developed for fits using censored data but the quality of fit can be judged using Akaike and Schwarz's Bayesian information criteria (AIC and BIC) and the goodness-of-fit CDF plot, respectively provided when summarizing or plotting an object of class 'fitdistcens'. Function `cdfcompcens` can also be used to compare the fit of various distributions to the same censored data set. Its call is similar to the one of `cdfcomp`. Below is an example of comparison of the two fitted distributions to the `salinity` data set (see Figure 10).

R> cdfcompcens(list(fsal.ln, fsal.l1), legendtext = c("lognormal",
+ "loglogistic"))
<figure><img src="image_13.png" /><figcaption>Figure 10: Goodness-of-fit CDF plot for fits of a lognormal and a loglogistic distribution to censored data: LC50 values from the salinity data set.</figcaption></figure>

Function `bootdistcens` is the equivalent of `bootdist` for censored data, except that it only proposes nonparametric bootstrap. Indeed, it is not obvious to simulate censoring within a parametric bootstrap resampling procedure. The generic function `quantile` can also be applied to an object of class 'fitdistcens' or 'bootdistcens', as for continuous non-censored data.

In addition to the fit of distributions to censored or non censored continuous data, our package can also accommodate discrete variables, such as count numbers, using the functions developed for continuous non-censored data. These functions will provide somewhat different graphs and statistics, taking into account the discrete nature of the modeled variable. The discrete nature of the variable is automatically recognized when a classical distribution is fitted to data (binomial, negative binomial, geometric, hypergeometric and Poisson distributions) but must be indicated by fixing argument `discrete` to `TRUE` in the call to functions in other cases. The `toxocara` data set included in the package corresponds to the observation of such a discrete variable. Numbers of *Toxocara cati* parasites present in the digestive tract are reported from a random sampling of feral cats living on Kerguelen island (Fromont, Morvilliers, Artois, and Pontier 2001). We will use it to illustrate the case of discrete data.

```R
R> data("toxocara", package = "fitdistrplus")
```

```R
R> str(toxocara)
```

```
/data.frame':          53 obs. of 1 variable:
```

```
$ number: int 0000000000 ...
```
<figure><img src="image_14.png" /><figcaption>Figure 11: Plot of the fit of a discrete distribution (Poisson distribution fitted to numbers of *Toxocara cati* parasites from the *toxocara* data set).</figcaption></figure>

The fit of a discrete distribution to discrete data by maximum likelihood estimation requires the same procedure as for continuous non-censored data. As an example, using the *toxocara* data set, Poisson and negative binomial distributions can be easily fitted.

R> (ftoxo.P <- fitdist(toxocara$number, "pois"))

Fitting of the distribution ' pois ' by maximum likelihood

Parameters:

<table><thead><tr><th></th><th>estimate</th><th>Std. Error</th></tr></thead><tbody><tr><td>lambda</td><td>8.679</td><td>0.4047</td></tr></tbody></table>

R> (ftoxo.nb <- fitdist(toxocara$number, "nbinom"))

Fitting of the distribution ' nbinom ' by maximum likelihood

Parameters:

<table><thead><tr><th></th><th>estimate</th><th>Std. Error</th></tr></thead><tbody><tr><td>size</td><td>0.3971</td><td>0.08289</td></tr><tr><td>mu</td><td>8.6803</td><td>1.93501</td></tr></tbody></table>

For discrete distributions, the plot of an object of class 'fitdist' simply provides two goodness-of-fit plots comparing empirical and theoretical distributions in density and in CDF (Figure 11). Function cdfcomp can also be used to compare several plots to the same data set, as follows for the previous fits (Figure 12).

R> plot(ftoxo.P)
R> cdfcomp(list(ftoxo.P, ftoxo.nb), legendtext = c("Poisson",
+ "negative binomial"))

When fitting discrete distributions, the Chi-squared statistic is computed by the `gofstat` function using cells defined by the argument `chisqbreaks` or cells automatically defined from
<figure><img src="image_15.png" /><figcaption>Figure 12: Comparison of the fits of a negative binomial and a Poisson distribution to numbers of *Toxocara cati* parasites from the *toxocara* data set.</figcaption></figure>

the data in order to reach roughly the same number of observations per cell. This number is roughly equal to the argument *meancount*, or slightly greater if there are some ties. The choice to define cells from the empirical distribution (data), and not from the theoretical distribution, was made to enable the comparison of Chi-squared values obtained with different distributions fitted on the same data set. If arguments *chisqbreaks* and *meancount* are both omitted, *meancount* is fixed in order to obtain roughly $(4n)^{2/5}$ cells, with *n* the length of the data set (Vose 2010). Using this default option the two previous fits are compared as follows, giving the preference to the negative binomial distribution, from both Chi-squared statistics and information criteria:

```
R> gofstat(list(ftoxo.P, ftoxo.nb),
+   fitnames = c("Poisson", "negative binomial"))
```

Chi-squared statistic: 31257 7.486
Degree of freedom of the Chi-squared distribution: 5 4
Chi-squared p-value: 0 0.1123
the p-value may be wrong with some theoretical counts < 5
Chi-squared table:

<table><thead><tr><th></th><th></th><th>obscounts</th><th>theo</th><th>Poisson</th><th>theo negative binomial</th></tr></thead><tbody><tr><td>&lt;= 0</td><td>14</td><td>0.009014</td><td></td><td>15.295</td><td></td></tr><tr><td>&lt;= 1</td><td>8</td><td>0.078237</td><td></td><td>5.809</td><td></td></tr><tr><td>&lt;= 3</td><td>6</td><td>1.321767</td><td></td><td>6.845</td><td></td></tr><tr><td>&lt;= 4</td><td>6</td><td>2.131298</td><td></td><td>2.408</td><td></td></tr><tr><td>&lt;= 9</td><td>6</td><td>29.827829</td><td></td><td>7.835</td><td></td></tr><tr><td>&lt;= 21</td><td>6</td><td>19.626224</td><td></td><td>8.271</td><td></td></tr><tr><td>&gt; 21</td><td>7</td><td>0.005631</td><td></td><td>6.537</td><td></td></tr></tbody></table>
Goodness-of-fit criteria

<table><thead><tr><th></th><th>Poisson negative binomial</th><th></th></tr></thead><tbody><tr><td>Aikake's Information Criterion</td><td>1017</td><td>322.7</td></tr><tr><td>Bayesian Information Criterion</td><td>1019</td><td>326.6</td></tr></tbody></table>

4. Conclusion

The R package fitdistrplus allows to easily fit distributions. Our main objective while developing this package was to provide tools for helping R users to fit distributions to data. We have been encouraged to pursue our work by feedbacks from users of our package in various areas as food or environmental risk assessment, epidemiology, ecology, molecular biology, genomics, bioinformatics, hydraulics, mechanics, financial and actuarial mathematics or operations research. Indeed, this package is already used by a lot of practitioners and academics for

* simple MLE fits (Jaloustre, Cornu, Morelli, Noel, and Delignette-Muller 2011; Sak and Haksoz 2011; Koch, Yemshanov, Magarey, and Smith 2012; Marquetoux et al. 2012; Scholl, Nice, Fordyce, Gompert, and Forister 2012; Suuronen, Kallonen, Eik, Puttonen, Serimaa, and Herrmann 2012; González-Varo, López-Bao, and Guitián 2012; Mandl, Monteiro, Vrisekoop, and Germain 2013; Malá 2013; Nadarajah and Bakar 2013; Simó, Casaña, and Sabaté 2013; Westphal-Fitch and Fitch 2013; Wayland 2013; Viana, Santamaría, Michot, and Figuerola 2013; Wu 2013a; Stagge and Moglen 2013; Fiorelli et al. 2013; Tikole, Jaravine, Orekhov, and Guentert 2013; Voigt et al. 2014),

* for MLE fits and goodness-of-fit statistics (Tarnczi, Fenyves, and Bcs 2011; Bagaria, Jaravine, Huang, Montelione, and Güntert 2012; Benavides-Piccione, Fernaud-Espinosa, Robles, Yuste, and DeFelipe 2012; Breitbach, Böhning-Gaese, Laube, and Schleuning 2012; Pouillot and Delignette-Muller 2010; Vaninsky 2013),

* for MLE fits and bootstrap (Croucher, Harris, Barquist, Parkhill, and Bentley 2012; Méheust, Cann, Reponen, Wakefield, and Vesper 2012; Orellano, Reynoso, Grassi, Palmieri, Uez, and Carlino 2012; Tello, Austin, and Telfer 2012; Hoelzer, Pouillot, Gallagher, Silverman, Kause, and Dennis 2012; Prosser, Hungerford, Erwin, Ottinger, Takekawa, and Ellis 2013; Zhang, Habib, Kuligowski, and Kim 2013; Rigaux, é, Albert, and Carlin 2014),

* for MLE fits, bootstrap and goodness-of-fit statistics (Larras, Montuelle, and Bouchez 2013),

* for MME fit (Luangkesorn, Norman, Zhuang, Falbo, and Sysko 2012; Callau Poduje, Belli, and Haberlandt 2013; Sato et al. 2013),

* for censored MLE and bootstrap (Leha et al. 2011; Pouillot, Hoelzer, Chen, and Dennis 2012; Jongenburger, Reij, Boer, Zwietering, and Gorris 2012; Commeau et al. 2012; De La Huerta Contreras, Vaquera Huerta, and Arnold 2015),

* for graphic analyzing in Anand, Yeturu, and Chandra (2012),
* for grouped-data fitting methods (Fu, Steiner, and Costafreda 2012) or

* more generally (Busschaert et al. 2010; Eling 2012; Samuel-Rosa, Dalmolin, and Miguel 2013; Srinivasan, Sorrell, Brooks, Edwards, and McDougle 2013; Meyer, Zhang, Hayakawa, Imai, and Przeworski 2013; Guillier, Danan, Bergis, Delignette-Muller, Granier, Rudelle, Beaufort, and Brisabois 2013; Daelman, Membré, Jacxsens, Vermeulen, Devlieghere, and Uyttendaele 2013; Eik, Luhmus, Tigasson, Listak, Puttonen, and Herrmann 2013; Wu 2013b; Drake, Chalabi, and Coker 2015).

The fitdistrplus package is complementary with the distrMod package (Kohl and Ruckdeschel 2010). distrMod provides an even more flexible way to estimate distribution parameters but its use requires a greater initial investment to learn how to manipulate the S4 classes and methods developed in the distr family packages. Many extensions of the fitdistrplus package are planned in the future: we target to extend to censored data some methods for the moment only available for non-censored data, especially concerning goodness-of-fit evaluation and fitting methods. We will also enlarge the choice of fitting methods for non-censored data, by proposing new goodness-of-fit distances (e.g., distances based on quantiles) for maximum goodness-of-fit estimation and new types of moments (e.g., limited expected values) for moment matching estimation. At last, we will consider the case of multivariate distribution fitting.

## Acknowledgments

The package would not have been at this stage without the stimulating contribution of Régis Pouillot and Jean-Baptiste Denis, especially for its conceptualization. We also want to thank Régis Pouillot for his very valuable comments on the first version of this paper. The authors gratefully acknowledges the two anonymous referees and the Editor for useful and constructive comments. The remaining errors, of course, should be attributed to the authors alone.

## References

* Anand P, Yeturu K, Chandra N (2012). “PocketAnnotate: Towards Site-Based Function Annotation.” *Nucleic Acids Research*, **40**(W1), W400–W408.

* Asquith WH (2014). *lmomco* – L-Moments, Trimmed L-Moments, L-Comoments, Censored L-Moments, and Many Distributions. R package version 2.1.1, URL http://CRAN.R-project.org/package=lmomco.

* Bagaria A, Jaravine V, Huang YJ, Montelione GT, Güntert P (2012). “Protein Structure Validation by Generalized Linear Model Root-Mean-Square Deviation Prediction.” *Protein Science*, **21**(2), 229–238.

* Benavides-Piccione R, Fernaud-Espinosa I, Robles V, Yuste R, DeFelipe J (2012). “Age-Based Comparison of Human Dendritic Spine Structure Using Complete Three-Dimensional Reconstructions.” *Cerebral Cortex*, **23**(8), 1798–1810.

* Blom G (1959). *Statistical Estimates and Transformed Beta Variables*. John Wiley & Sons.
Breitbach N, Böhning-Gaese K, Laube I, Schleuning M (2012). “Short Seed-Dispersal Distances and Low Seedling Recruitment in Farmland Populations of Bird-Dispersed Cherry Trees.” *Journal of Ecology*, **100**(6), 1349–1358.

Busschaert P, Geeraerd AH, Uyttendaele M, VanImpe JF (2010). “Estimating Distributions out of Qualitative and (Semi)Quantitative Microbiological Contamination Data for Use in Risk Assessment.” *International Journal of Food Microbiology*, **138**(3), 260–269.

Callau Poduje AC, Belli A, Haberlandt U (2013). “Dam Risk Assessment Based on Univariate versus Bivariate Statistical Approaches: A Case Study for Argentina.” *Hydrological Sciences Journal*, **59**(12), 2216–2232.

Casella G, Berger RL (2002). *Statistical Inference*. 2nd edition. Duxbury Thomson Learning.

Commeau N, Parent E, Delignette-Muller ML, Cornu M (2012). “Fitting a Lognormal Distribution to Enumeration and Absence/Presence Data.” *International Journal of Food Microbiology*, **155**(3), 146–152.

Croucher NJ, Harris SR, Barquist L, Parkhill J, Bentley SD (2012). “A High-Resolution View of Genome-Wide Pneumococcal Transformation.” *PLoS Pathogens*, **8**(6), e1002745.

Cullen AC, Frey HC (1999). *Probabilistic Techniques in Exposure Assessment*. 1st edition. Plenum Publishing Co.

Daelman J, Membré JM, Jacxsens L, Vermeulen A, Devlieghere F, Uyttendaele M (2013). “A Quantitative Microbiological Exposure Assessment Model for *Bacillus Cereus* in REPFEDs.” *International Journal of Food Microbiology*, **166**(3), 433–449.

D'Agostino RB, Stephens MA (1986). *Goodness-of-Fit Techniques*. 1st edition. Dekker.

De La Huerta Contreras V, Vaquera Huerta H, Arnold BC (2015). “A Test for Equality of Variance with Censored Samples.” *Journal of Statistical Computation and Simulation*, **85**(3), 450–467.

Delignette-Muller ML, Cornu M, AFSSA-STEC-Study-Group (2008). “Quantitative Risk Assessment for *Escherichia Coli* O157:H7 in Frozen Ground Beef Patties Consumed by Young Children in French Households.” *International Journal of Food Microbiology*, **128**(1), 158–164.

Delignette-Muller ML, Pouillot R, Denis JB, Dutang C (2015). ***fitdistrplus: Help to Fit of a***
***Parametric Distribution to Non-Censored or Censored Data.*** R package version 1.0-4, URL http://CRAN.R-project.org/package=fitdistrplus.

Drake T, Chalabi Z, Coker R (2015). “Buy Now, Saved Later? The Critical Impact of Time-to-Pandemic Uncertainty on Pandemic Cost-Effectiveness Analyses.” *Health Policy and Planning*, **30**(1), 100–110.

Dutang C (2014). “CRAN Task View: Probability Distributions.” Version 2014-12-18, URL http://CRAN.R-project.org/view=Distributions.

Dutang C, Goulet V, Pigeon M (2008). “**actuar**: An R Package for Actuarial Science.” *Journal of Statistical Software*, **25**(7), 1–37. URL http://www.jstatsoft.org/v25/i07/.
Efron B, Tibshirani RJ (1994). *An Introduction to the Bootstrap*. 1st edition. Chapman & Hall.

Eik M, Luhmus K, Tigasson M, Listak M, Puttonen J, Herrmann H (2013). “DC-Conductivity Testing Combined with Photometry for Measuring Fibre Orientations in SFRC.” *Journal of Materials Science*, **48**(10), 3745–3759.

Eling M (2012). “Fitting Insurance Claims to Skewed Distributions: Are the Skew-Normal and the Skew-Student Good Models?” *Insurance: Mathematics and Economics*, **51**(2), 239–248.

Fiorelli LE, Ezcurra MD, Hechenleitner EM, Argañaraz E, Jeremias R, Taborda A, Trotteyn MJ, von Baczko MB, Desojo JB (2013). “The Oldest Known Communal Latrines Provide Evidence of Gregarism in Triassic Megaherbivores.” *Scientific Reports*, **3**(3348), 1–7.

Fromont E, Morvilliers L, Artois M, Pontier D (2001). “Parasite Richness and Abundance in Insular and Mainland Feral Cats: Insularity or Density?” *Parasitology*, **123*(Part 2), 143–151.

Fu CHY, Steiner H, Costafreda SG (2012). “Predictive Neural Biomarkers of Clinical Response in Depression: A Meta-Analysis of Functional and Structural Neuroimaging Studies of Pharmacological and Psychological Therapies.” *Neurobiology of Disease*, **52**, 75–83.

González-Varo JP, López-Bao JV, Guitián J (2012). “Functional Diversity Among Seed Dispersal Kernels Generated by Carnivorous Mammals.” *Journal of Animal Ecology*, **82**, 562–571.

Goulet V (2014). *actuar: An R Package for Actuarial Science*. R package version 1.1-8, URL http://CRAN.R-project.org/package=actuar.

Guillier L, Danan C, Bergis H, Delignette-Muller ML, Granier S, Rudelle S, Beaufort A, Brisabois A (2013). “Use of Quantitative Microbial Risk Assessment when Investigating Foodborne Illness Outbreaks: The Example of a Monophasic *Salmonella Typhimurium* 4,5,12:i- Outbreak Implicating Beef Burgers.” *International Journal of Food Microbiology*, **166**(3), 471–478.

Helsel DR (2005). *Nondetects and Data Analysis: Statistics for Censored Environmental Data*. John Wiley & Sons.

Hirano SS, Clayton MK, Upper CD (1994). “Estimation of and Temporal Changes in Means and Variances of Populations of *Pseudomonas Syringae* on Snap Bean Leaflets.” *Phytopathology*, **84**(9), 934–940.

Hoelzer K, Pouillot R, Gallagher D, Silverman MB, Kause J, Dennis S (2012). “Estimation of *Listeria Monocytogenes* Transfer Coefficients and Efficacy of Bacterial Removal through Cleaning and Sanitation.” *International Journal of Food Microbiology*, **157**(2), 267–277.

Hose GC, Van den Brink PJ (2004). “Confirming the Species-Sensitivity Distribution Concept for Endosulfan Using Laboratory, Mesocosm, and Field Data.” *Archives of Environmental Contamination and Toxicology*, **47**(4), 511–520.
Hyndman RJ, Fan Y (1996). "Sample Quantiles in Statistical Packages." *The American Statistician*, **50**(4), 361-365.

Jaloustre S, Cornu M, Morelli E, Noel V, Delignette-Muller ML (2011). "Bayesian Modeling of *Clostridium Perfringens* Growth in Beef-in-Sauce Products." *Food Microbiology*, **28**(2), 311-320.

Jongenburger I, Reij MW, Boer EPJ, Zwietering MH, Gorris LGM (2012). "Modelling Heterogeneous and Heterogeneous Microbial Contaminations in a Powdered Food Product." *International Journal of Food Microbiology*, **157**(1), 35-44.

Jordan D (2005). "Simulating the Sensitivity of Pooled-Sample Herd Tests for Fecal Salmonella in Cattle." *Preventive Veterinary Medicine*, **70**(1-2), 59-73.

Karvanen J (2006). "Estimation of Quantile Mixtures via L-Moments and Trimmed L-Moments." *Computational Statistics & Data Analysis*, **51**(2), 947-959.

Karvanen J (2014). *Lmoments: L-Moments and Quantile Mixtures*. R package version 1.1-6, URL http://CRAN.R-project.org/package=Lmoments.

Kefford BJ, Fields EJ, Clay C, Nugegoda D (2007). "Salinity Tolerance of Riverine Macroinvertebrates from the Southern Murray-Darling Basin." *Marine and Freshwater Research*, **58**(11), 1019-1031.

Klein JP, Moeschberger ML (2003). *Survival Analysis: Techniques for Censored and Truncated Data*. 2nd edition. Springer-Verlag.

Klugman SA, Panjer HH, Willmot GE (2009). *Loss Models: From Data to Decisions*. 3rd edition. John Wiley & Sons.

Koch FH, Yemshanov D, Magarey RD, Smith WD (2012). "Dispersal of Invasive Forest Insects via Recreational Firewood: A Quantitative Analysis." *Journal of Economic Entomology*, **105**(2), 438-450.

Kohl M, Ruckdeschel P (2010). "R Package distrMod: S4 Classes and Methods for Probability Models." *Journal of Statistical Software*, **35**(10), 1-27. URL http://www.jstatsoft.org/v35/i10/.

Larras F, Montuelle B, Bouchez A (2013). "Assessment of Toxicity Thresholds in Aquatic Environments: Does Benthic Growth of Diatoms Affect Their Exposure and Sensitivity to Herbicides?" *Science of the Total Environment*, **463-464**, 469-477.

Leha A, Beissbarth T, Jung K (2011). "Sequential Interim Analyses of Survival Data in DNA Microarray Experiments." *BMC Bioinformatics*, **12**(127), 1-14.

Luangkesorn KL, Norman BA, Zhuang Y, Falbo M, Sysko J (2012). "Practice Summaries: Designing Disease Prevention and Screening Centers in Abu Dhabi." *Interfaces*, **42**(4), 406-409.

Luceno A (2006). "Fitting the Generalized Pareto Distribution to Data Using Maximum Goodness-of-Fit Estimators." *Computational Statistics & Data Analysis*, **51**(2), 904-917.
Malá I (2013). "The Use of Finite Mixtures of Lognormal and Gamma Distributions." *Research Journal of Economics, Business and ICT*, **8**(2), 55–61.

Mandl JN, Monteiro JP, Vrisekoop N, Germain RN (2013). "T Cell-Positive Selection Uses Self-Ligand Binding Strength to Optimize Repertoire Recognition of Foreign Antigens." *Immunity*, **38**(2), 263–274.

Marquetoux N, Paul M, Wongnarkpet S, Poolkhet C, Thanapongtham W, Roger F, Ducrot C, Chalvet-Monfray K (2012). "Estimating Spatial and Temporal Variations of the Reproduction Number for Highly Pathogenic Avian Influenza H5N1 Epidemic in Thailand." *Preventive Veterinary Medicine*, **106**(2), 143–151.

McNeil AJ (1997). "Estimating the Tails of Loss Severity Distributions Using Extreme Value Theory." *ASTIN Bulletin*, **27**(1), 117–137.

Mebane, Jr WR, Sekhon JS (2011). "Genetic Optimization Using Derivatives: The rgenoud Package for R." *Journal of Statistical Software*, **42**(11), 1–26. URL http://www.jstatsoft.org/v42/i11/.

Méheust D, Cann PL, Reponen T, Wakefield J, Vesper S (2012). "Possible Application of the Environmental Relative Moldiness Index in France: a Pilot Study in Brittany." *International Journal of Hygiene and Environmental Health*, **216**(3), 333–340.

Meyer WK, Zhang S, Hayakawa S, Imai H, Przeworski M (2013). "The Convergent Evolution of Blue Iris Pigmentation in Primates Took Distinct Molecular Paths." *American Journal of Physical Anthropology*, **151**(3), 398–407.

Nadarajah S, Bakar SAA (2013). "CompLognormal: An R Package for Composite Lognormal Distributions." *The R Journal*, **5**(2), 98–104.

Orellano PW, Reynoso JI, Grassi A, Palmieri A, Uez O, Carlino O (2012). "Estimation of the Serial Interval for Pandemic Influenza (pH1N1) in the Most Southern Province of Argentina." *Iranian Journal of Public Health*, **41**(12), 26–29.

Poncet P (2012). *modeest: Mode Estimation*. R package version 2.1, URL http://CRAN.R-project.org/package=modeest.

Posthuma L, Suter GW, Traas TP (2010). *Species Sensitivity Distributions in Ecotoxicology*. Environmental and Ecological Risk Assessment Series. Taylor & Francis.

Pouillot R, Delignette-Muller ML (2010). "Evaluating Variability and Uncertainty Separately in Microbial Quantitative Risk Assessment Using two R Packages." *International Journal of Food Microbiology*, **142**(3), 330–340.

Pouillot R, Delignette-Muller ML, Denis JB (2011). *mc2d: Tools for Two-Dimensional Monte-Carlo Simulations*. R package version 0.1-12, URL http://CRAN.R-project.org/package=mc2d.

Pouillot R, Hoelzer K, Chen Y, Dennis S (2012). "Estimating Probability Distributions of Bacterial Concentrations in Food Based on Data Generated Using the Most Probable Number (MPN) Method for Use in Risk Assessment." *Food Control*, **29**(2), 350–357.
Prosser DJ, Hungerford LL, Erwin RM, Ottinger MA, Takekawa JY, Ellis EC (2013). “Mapping Avian Influenza Transmission Risk at the Interface of Domestic Poultry and Wild Birds.” *Frontiers in Public Health*, *1*(28), 1–11.

R Core Team (2014). *R: A Language and Environment for Statistical Computing*. R Foundation for Statistical Computing, Vienna, Austria. URL http://www.R-project.org/.

Ricci V (2005). “Fitting Distributions with R.” Contributed documentation available on CRAN, URL http://CRAN.R-project.org/doc/contrib/Ricci-distributions-en.pdf.

Rigaux C, é SA, Albert I, Carlin F (2014). “Quantitative Assessment of the Risk of Microbial Spoilage in Foods. Prediction of Non-Stability at 55°C Caused by *Geobacillus Stearothermophilus* in Canned Green Beans.” *International Journal of Food Microbiology*, *171*, 119–128.

Sak H, Haksoz C (2011). “A Copula-Based Simulation Model for Supply Portfolio Risk.” *Journal of Operational Risk*, *6*(3), 15–38.

Samuel-Rosa A, Dalmolin RSD, Miguel P (2013). “Building Predictive Models of Soil Particle-Size Distribution.” *Revista Brasileira de Ciencia do Solo*, **37**(2), 422–430.

Sato MIZ, Galvani AT, Padula JA, Nardocci AC, de Souza Lauretto M, Razzolini MTP, Hachich EM (2013). “Assessing the Infection Risk of *Giardia* and *Cryptosporidium* in Public Drinking Water Delivered by Surface Water Systems in Sao Paulo State, Brazil.” *Science of the Total Environment*, **442**, 389–396.

Scholl CF, Nice CC, Fordyce JA, Gompert Z, Forister ML (2012). “Larval Performance in the Context of Ecological Diversification and Speciation in Lycaeides Butterflies.” *International Journal of Ecology*, **2012**(ID 242154), 1–13.

Simó J, Casaña F, Sabaté J (2013). “Modelling “Calçots” (*Alium Cepa L.*) Growth by Gompertz Function.” *Statistics and Operations Research Transactions*, **37**(1), 95–106.

Srinivasan S, Sorrell TP, Brooks JP, Edwards DJ, McDougle RD (2013). “Workforce Assessment Method for an Urban Police Department: Using Analytics to Estimate Patrol Staffing.” *Policing: An International Journal of Police Strategies & Management*, **36**(4), 702–718.

Stagge JH, Moglen GE (2013). “A Nonparametric Stochastic Method for Generating Daily Climate-Adjusted Streamflows.” *Water Resources Research*, **49**(10), 6179–6193.

Suuronen JP, Kallonen A, Eik M, Puttonen J, Serimaa R, Herrmann H (2012). “Analysis of Short Fibres Orientation in Steel Fibre-Reinforced Concrete (SFRC) by X-Ray Tomography.” *Journal of Materials Science*, **48**(3), 1358–1367.

Tarnczi T, Fenyves V, Bcs Z (2011). “The Business Uncertainty and Variability Management with Real Options Models Combined two Dimensional Simulation.” *International Journal of Management Cases*, **13**(3), 159–167.

Tello A, Austin B, Telfer TC (2012). “Selective Pressure of Antibiotic Pollution on Bacteria of Importance to Public Health.” *Environmental Health Perspectives*, **120**(8), 1100–1106.
Therneau TM (2014). *survival: Survival Analysis*. R package version 2.37-7, URL http://CRAN.R-project.org/package=survival.

Therneau TM, Grambsch PM (2000). *Modeling Survival Data: Extending the Cox Model*. Springer-Verlag, New York.

Tikole S, Jaravine V, Orekhov VY, Guentert P (2013). "Effects of NMR Spectral Resolution on Protein Structure Calculation." *PloS ONE*, **8**(7), e68567.

Tse YK (2009). *Nonlife Actuarial Models: Theory, Methods and Evaluation*. International Series on Actuarial Science, 1st edition. Cambridge University Press.

Turnbull BW (1974). "Nonparametric Estimation of a Survivorship Function with Doubly Censored Data." *Journal of the American Statistical Association*, **69**(345), 169-173.

Vaninsky AY (2013). "Stochastic DEA with a Perfect Object and Its Application to Analysis of Environmental Efficiency." *American Journal of Applied Mathematics and Statistics*, **1**(4), 57-63.

Venables WN, Ripley BD (2010). *Modern Applied Statistics with S*. 4th edition. Springer-Verlag.

Viana DS, Santamaría L, Michot TC, Figuerola J (2013). "Allometric Scaling of Long-Distance Seed Dispersal by Migratory Birds." *The American Naturalist*, **181**(5), 649-662.

Voigt CC, Lehnert LS, Popa-Lisseanu AG, Ciechanowski M, Estók P, Gloza-Rausch F, Goerfoel T, Goettsche M, Harrje C, Hoetzel M, Teige T, Wohlgemuth R, Kramer-Schadt S (2014). "The Trans-Boundary Importance of Artificial Bat *Hibernacula* in Managed European Forests." *Biodiversity and Conservation*, **23**, 617-631.

Vose D (2010). *Quantitative Risk Analysis. A Guide to Monte Carlo Simulation Modelling*. 1st edition. John Wiley & Sons.

Wayland MT (2013). "Morphological Variation in *Echinorhynchus truttae* Schrank, 1788 and the *Echinorhynchus bothniensis* Zdzitowiecki & Valtonen, 1987 Species Complex from Freshwater Fishes of Northern Europe." *Biodiversity Data Journal*, **1**, e975.

Westphal-Fitch G, Fitch WT (2013). "Spatial Analysis of "Crazy Quilts", A Class of Potentially Random Aesthetic Artefacts." *PloS ONE*, **8**(9), e74055.

Wu XZ (2013a). "Probabilistic Slope Stability Analysis by a Copula-Based Sampling Method." *Computational Geosciences*, **17**(5), 739-755.

Wu XZ (2013b). "Trivariate Analysis of Soil Ranking-Correlated Characteristics and Its Application to Probabilistic Stability Assessments in Geotechnical Engineering Problems." *Soils and Foundations*, **53**(4), 540-556.

Zhang Y, Habib E, Kuligowski RJ, Kim D (2013). "Joint Distribution of Multiplicative Errors in Radar and Satellite QPEs and its Use in Estimating the Conditional Exceedance Probability." *Advances in Water Resources*, **59**, 133-145.
Affiliation:

Marie Laure Delignette-Muller

Université de Lyon

Université Lyon 1, CNRS, UMR5558, Laboratoire de Biométrie et Biologie évolutive

VetAgro Sup, Campus Vétérinaire de Lyon

1, avenue Bourgelat

69820 Marcy L'Etoile, France

E-mail: marielaure.delignettemuller@vetagro-sup.fr

URL: http://lbbe.univ-lyon1.fr/~Delignette-Muller-Marie-Laure-.html