This function fits 43 different continuous distributions to the Monte Carlo simulation results obtained by propagate, or to any other vector of observations, and ranks them by ascending BIC. By default, the distributions are fitted by maximum likelihood to the observations; the (weighted) least-squares histogram fit used up to version 1.2-0 is available with method = "nls".
fitDistr(object, method = c("mle", "nls"), nbin = 100, nmax = 1E5,
weights = FALSE, verbose = TRUE, brute = c("fast", "slow"),
plot = c("hist", "qq", "none"), distsel = NULL, ...)A list with the following items:
stat: a data frame of the fitted distributions, sorted by ascending BIC. For method = "mle": the number of parameters k, logLik, AIC, BIC, dBIC (difference to the lowest BIC) and the Schwarz weight; distributions that could not be fitted, or whose support excludes the data, have NA. For method = "nls": BIC, the (weighted) residual sum-of-squares RSS and the mean squared residual MSE.
fit: a list with the fit for each distribution, sorted as in stat. For "mle": a list with par, se, logLik, k, n, convergence (of optim) and the internal parameters theta; for "nls": the minpack.lm::nls.lm object.
par: a list of the estimated parameters, named as the arguments of the density functions.
se: a list of the parameters' standard errors (NA for non-regular parameters with "mle").
dens: a list with the density functions of the fitted distributions, sorted and named as in fit.
bestfit: the best model in terms of lowest BIC.
bestpar: the parameters of bestfit.
bestse: the parameters' standard errors of bestfit.
fitted: the density of bestfit at the histogram bin midpoints.
residuals: the histogram density minus fitted.
method: the fitting method.
n: the number of observations used (for "mle", at most nmax).
an object of class 'propagate' or a vector containing observations.
"mle" (default) for maximum likelihood estimation on the observations, "nls" for the least-squares fit of the density to the histogram. See 'Details'.
the number of bins in the histogram. For method = "mle", the histogram is only used for the plots and for fitted/residuals.
method = "mle" only: the maximum number of observations used for fitting. Larger data sets, such as propagate output with \(10^6\) values, are reduced to an evenly spaced, deterministic subsample of nmax values.
method = "nls" only: numeric or logical. Either a numeric vector of weights with one weight per histogram bin, or if TRUE, the distributions are fitted with weights = 1/(counts per bin); bins with zero counts get weight 1. Ignored, with a message, for method = "mle".
logical. If TRUE, steps of the analysis are printed to the console.
complexity of the brute force approach. See 'Details'.
if "hist", a plot with the "best" distribution (in terms of lowest BIC) on top of the histogram is displayed. If "qq", the fitted density values of the best distribution are plotted against the histogram densities, bin by bin, with the identity line (a density-density plot, not a quantile-quantile plot).
a vector of distribution numbers to select from the complete cohort as listed below, e.g. c(1:10, 15). All outputs (fit, par, se, dens) refer to the selected distributions, sorted by BIC.
other parameters to be passed to the plots.
Andrej-Nikolai Spiess
Fits the following 43 distributions (definitions in Johnson, Kotz & Balakrishnan 1994, 1995 and Leemis & McQueston 2008, unless another reference is given; parametrisations as in the density functions of the package):
1) Normal distribution (dnorm)
2) Skewed-normal distribution (propagate:::dsn)
3) Generalized normal distribution (propagate:::dgnorm)
4) Log-normal distribution (dlnorm)
5) Scaled and shifted t-distribution (propagate:::dst) => JCGM 101:2008, 6.4.9.2. The parameter sd is the scale, not the standard deviation.
6) Logistic distribution (dlogis)
7) Uniform distribution (dunif)
8) Triangular distribution (propagate:::dtriang)
9) Trapezoidal distribution (propagate:::dtrap)
10) Curvilinear Trapezoidal distribution (propagate:::dctrap) => JCGM 101:2008, 6.4.3
11) Gamma distribution (dgamma)
12) Inverse Gamma distribution (propagate:::dinvgamma)
13) Cauchy distribution (dcauchy)
14) Laplace distribution with standard deviation sigma (propagate:::dlaplace)
15) Gumbel distribution (propagate:::dgumbel)
16) Johnson SU distribution (propagate:::dJSU)
17) Johnson SB distribution (propagate:::dJSB)
18) Three-parameter Weibull distribution (propagate:::dweibull2)
19) Two-parameter beta distribution (dbeta)
20) Four-parameter beta distribution (propagate:::dbeta2)
21) Arcsine distribution (propagate:::darcsin)
22) Von Mises distribution (propagate:::dmises)
23) Inverse Gaussian distribution (propagate:::dinvgauss)
24) Generalized Extreme Value distribution (propagate:::dgevd)
25) Rayleigh distribution (propagate:::drayleigh)
26) Chi-square distribution (dchisq)
27) Exponential distribution (dexp)
28) F-distribution (df)
29) Burr type VIII distribution (propagate:::dburr) => Johnson, Kotz & Balakrishnan (1994)
30) Chi distribution (propagate:::dchi)
31) Inverse Chi-square distribution (propagate:::dinvchisq)
32) Raised cosine distribution (propagate:::dcosine)
33) Pareto distribution (propagate:::dpareto)
34) Levy distribution (propagate:::dlevy)
35) Gompertz distribution (propagate:::dgompertz)
36) 3P log-normal distribution, Johnson SL (propagate:::dlnorm3): \(\xi + \exp(\mu + \sigma Z)\), positive skewness with a free lower limit \(\xi\)
37) 3P gamma distribution, Pearson type III (propagate:::dgamma3): gamma distribution shifted by \(\xi\)
38) Rectangular-normal convolution (propagate:::dunorm): \(g(x) = [\Phi((x - a)/\sigma) - \Phi((x - b)/\sigma)]/(b - a)\), the output of models dominated by one rectangular input
39) Exponential power distribution (propagate:::dexpow): \(g(x) \propto \exp(-(|x - \mu|/\alpha)^\beta)\); Laplace for \(\beta = 1\), normal for \(\beta = 2\), rectangular for \(\beta \to \infty\)
40) Skew-t distribution (propagate:::dskewt, Azzalini & Capitanio 2003): skewness and heavy tails, e.g. for ratios of t-distributed inputs
41) Rice distribution (propagate:::drice): magnitude of a complex quantity with Gaussian real and imaginary parts; for large arguments, the exponentially scaled Bessel function \(I_0\) is replaced by its asymptotic expansion
42) Folded normal distribution (propagate:::dfnorm): absolute value of a normal variable
43) Two-component normal mixture (propagate:::dnmix): bimodal output distributions
All density functions return 0 outside their support (and -Inf for log = TRUE). Distributions whose support cannot contain the data (e.g. the log-normal for negative values, or the 2P beta outside \((0, 1)\)) are skipped. The density functions propagate:::d... are internal functions (hence :::). The code for the density functions, and the parametrisations used for maximum likelihood fitting (internal function .mleSpec), are in file "distr-densities.R".
Maximum likelihood (method = "mle", default).
For each distribution with density \(g(x \mid \theta)\), the parameters are estimated by maximising the log-likelihood of the \(n\) observations,
$$\ln L(\theta) = \sum_{i=1}^{n} \ln g(x_i \mid \theta),$$
and the distributions are ranked by
$$\mathrm{BIC} = -2\ln L(\hat\theta) + k \ln n, \qquad \mathrm{AIC} = -2\ln L(\hat\theta) + 2k,$$
with \(k\) the number of distribution parameters. Unlike the histogram fit, the BIC is then a proper likelihood-based criterion that does not depend on nbin. stat also reports \(\Delta\mathrm{BIC}\) to the best distribution and the Schwarz weights \(w_j = \exp(-\Delta\mathrm{BIC}_j/2)/\sum_l \exp(-\Delta\mathrm{BIC}_l/2)\), which approximate the posterior probabilities of the candidates: a weight close to 1 indicates a clear winner, several similar weights indicate that the data cannot distinguish between the candidates (e.g. nested families such as the normal and the t-distribution with large \(\nu\)). With large \(n\), the BIC discriminates strongly, and the best distribution is the best approximation among the candidates, which need not contain the true distribution.
Parametrisation. The likelihood is maximised over unconstrained internal parameters: scale and shape parameters on the log scale, modes and other interior proportions on the logit scale, and support limits as \(\mathrm{min}(x) - s - e^{\vartheta}\) and \(\mathrm{max}(x) + s + e^{\vartheta}\), so that the support always contains the data. Here \(s = (\max(x) - \min(x))/n\) is the mean spacing of the observations: keeping support limits at least one spacing outside the data bounds the likelihood of densities that are infinite at a support limit (arcsine, 3P Weibull with shape < 1), whose maximum likelihood would otherwise be infinite. For the uniform distribution, this gives practically the minimum-variance unbiased estimates \(\min(x) - (\max(x) - \min(x))/(n - 1)\) and \(\max(x) + (\max(x) - \min(x))/(n - 1)\) of the limits. Dependent parameters are mapped internally, e.g. the mode of the triangular distribution as \(b = a + (c - a)\,p\) with \(0 < p < 1\), or the upper limit \(\xi + \lambda\) of the Johnson SB distribution.
Optimisation. A grid of starting values around moment-based estimates is evaluated (brute = "fast": 3 values per parameter, \(3^k\) points; brute = "slow": 5 values, \(5^k\) points), with steps of a factor 10 for scale/shape parameters and of one standard deviation for location parameters. The likelihood is maximised by Nelder-Mead from the best grid point ("fast") or the best three ("slow"), followed by a BFGS refinement.
For the normal mixture 43), the component standard deviations are bounded below by 1 % of the sample standard deviation, because the likelihood of a mixture is unbounded when one component collapses onto a single value.
Standard errors. The standard errors are obtained from the inverse Hessian of \(-\ln L\) with respect to the internal parameters, transformed to the parameters of the density by the delta method. Parameters that are limits of the support are non-regular: their estimates converge to the sample extremes at rate \(1/n\) rather than \(1/\sqrt{n}\), and the Hessian does not describe their uncertainty. They are held fixed for the Hessian, and their standard errors are NA. This applies to the limits of the uniform, triangular, trapezoidal, curvilinear trapezoidal, 4P beta, arcsine and Johnson SB distributions, to the mode of the triangular and the inner corners of the trapezoidal distribution (kinks of the density), the location of the 3P Weibull and Rayleigh distributions, the scale of the Pareto distribution and \(\sigma\) of the cosine distribution.
Sample size. The likelihood is evaluated on all observations, or on an evenly spaced, deterministic subsample of nmax values (no random numbers are drawn). \(10^5\) values are ample to identify a distribution, and the run time is then similar to that of the histogram fit.
Histogram fit (method = "nls", versions up to 1.2-0).
The density is fitted by (weighted) least squares to the histogram densities at the bin midpoints, with minpack.lm::nls.lm. The parameter space is searched with a brute force approach, extended over three orders of magnitude \((0.1, 1, 10)\times \beta_i\) when brute = "fast", or five orders \((0.01, 0.1, 1, 10, 100)\times \beta_i\) when brute = "slow". The goodness-of-fit is calculated with BIC from the (weighted) Gaussian log-likelihood of the NLS fit, as in stats:::logLik.nls:
$$\ln(L) = \frac{1}{2}\left(\sum_{i=1}^{N} \ln(w_i) - N\left(\ln(2\pi) + 1 - \ln(N) + \ln\sum_{i=1}^{N} r_i^2\right)\right)$$
$$\mathrm{BIC} = -2\ln(L) + k \ln(N)$$
with \(r_i = \sqrt{w_i}\,(\hat{g}_i - g(x_i))\) the weighted residuals between the histogram density \(\hat{g}_i\) and the fitted density at the bin midpoints \(x_i\), \(w_i\) the weights, \(N\) the number of bins with non-zero weight, and \(k\) the number of distribution parameters plus one (for the residual variance). Because \(N\) is the number of bins, not of observations, these BIC values depend on nbin and are a heuristic criterion for comparing distributions on the same histogram. Setting weights = TRUE gives more weight to bins with low counts; also check whether the ranking remains stable when nbin is changed.
Use. Besides characterising the output distribution of propagate, fitDistr can identify a distribution for raw data prior to Monte Carlo simulation. A decent number of observations is needed to distinguish between similar distributions. In contrast to software that ranks by residual sum-of-squares or Kolmogorov-Smirnov/Anderson-Darling statistics, the BIC penalises the number of parameters and thus compensates for overfitting; the fits are maximum likelihood fits as in MASS::fitdistr.
Evaluation of measurement data - Supplement 1 to the "Guide to the expression of uncertainty in measurement" - Propagation of distributions using a Monte Carlo method.
JCGM 101:2008. Joint Committee for Guides in Metrology (2008).
tools:::Rd_expr_doi("10.59161/JCGM101-2008")
Distributions generated by perturbation of symmetry with emphasis on a multivariate skew t-distribution.
Azzalini A, Capitanio A.
Journal of the Royal Statistical Society B (2003), 65: 367-389.
tools:::Rd_expr_doi("10.1111/1467-9868.00391")
Continuous Univariate Distributions, Volume 1, 2nd edition.
Johnson NL, Kotz S, Balakrishnan N.
Wiley, New York (1994). ISBN 0-471-58495-9.
Continuous Univariate Distributions, Volume 2, 2nd edition.
Johnson NL, Kotz S, Balakrishnan N.
Wiley, New York (1995). ISBN 0-471-58494-0.
Univariate distribution relationships.
Leemis LM, McQueston JT.
The American Statistician (2008), 62: 45-53.
tools:::Rd_expr_doi("10.1198/000313008X270448")
Estimating the dimension of a model.
Schwarz G.
The Annals of Statistics (1978), 6: 461-464.
tools:::Rd_expr_doi("10.1214/aos/1176344136")
Model Selection and Multimodel Inference.
Burnham KP, Anderson DR.
Springer, New York, 2nd edition (2002).
tools:::Rd_expr_doi("10.1007/b97636")
## Linear example, small error
## => Normal distribution.
## Don't set 'distsel' for full test
EXPR1 <- expression(x + y)
x <- c(5, 0.1)
y <- c(1, 0.01)
DF1 <- cbind(x, y)
RES1 <- propagate(expr = EXPR1, data = DF1, type = "stat",
nsim = 10000, check = FALSE)
FIT1 <- fitDistr(RES1, distsel = 1:10)
FIT1$stat
# \donttest{
## Identify the distribution of raw data: t-distributed sample.
set.seed(1)
X <- propagate:::rst(10000, mean = 3, sd = 2, df = 5) # internal generator, see ?rDistr
FIT2 <- fitDistr(X, distsel = 1:10)
FIT2$stat # the t-distribution wins, weight close to 1
FIT2$bestpar; FIT2$bestse
# }
Run the code above in your browser using DataLab