propagate (version 1.3-0)

corDist: Correlated Monte Carlo samples of input quantities with arbitrary distributions

Description

Creates Monte Carlo samples of correlated input quantities, each with its own distribution, joined by a Gaussian copula: Gaussian, rectangular, triangular, arcsine, curvilinear trapezoidal or scaled and shifted t given by estimate and standard uncertainty, or any distribution given by a sampler (e.g. the generators of rDistr or of R's distributions) or by a sample of values. The result is used as input of propagate with type = "sim" and cov = "full".

Usage

corDist(mean, u, dist = "normal", cor = NULL, df = Inf, d = NA, nsim = 1E6,
        cor.type = c("pearson", "copula"))

Value

A matrix with nsim rows and one named column per input quantity, with the attributes "budget" (data frame: name, mean, u, df, dist) and "corMat.copula" (the copula correlation matrix \(\mathbf{R}_Z\) used).

Arguments

mean

a named vector of the estimates \(x_i\) of the input quantities; the names are the variables of the expression passed to propagate. NA for inputs given by a sampler or a sample.

u

the standard uncertainties \(u_i\), recycled to the length of mean (NA for samplers and samples).

dist

the distributions, a character vector or list (recycled): "normal", "rect", "tri" (symmetric triangular), "arcsine" (U-shaped), "ctrap" (curvilinear trapezoidal, JCGM 101:2008, 6.4.3, with d), "t" (scaled and shifted t-distribution with scale \(u_i\) and df degrees of freedom, JCGM 101:2008, 6.4.9), a sampler function(n) returning n random values, e.g. function(n) propagate:::rtrap(n, 0, 1, 2, 3) (the generators of rDistr are internal functions, hence :::) or function(n) rgamma(n, 2, 1), or a numeric vector of sampled or measured values (resampled with replacement if its length differs from nsim).

cor

the correlation matrix \(\mathbf{R}\) of the input quantities, in the order of mean; NULL for independent inputs.

df

the degrees of freedom of "t" margins (recycled). They are also passed to the GUM uncertainty framework of propagate.

d

the inexactness of the limits of "ctrap" margins (recycled).

nsim

the number of Monte Carlo samples.

cor.type

"pearson" (default): the copula correlation is chosen such that the Pearson correlation of the samples equals cor; "copula": cor is used as the parameter of the Gaussian copula. See 'Details'.

Author

Andrej-Nikolai Spiess

Details

Margins. Each input quantity has expectation \(x_i\) and standard deviation \(u_i\), with the half-widths \(a = \sqrt{3}\,u_i\) (rectangular), \(\sqrt{6}\,u_i\) (triangular), \(\sqrt{2}\,u_i\) (arcsine) and \(u_i^2 = a^2/3 + d^2/9\) (curvilinear trapezoidal). The t-margin has scale \(u_i\) and variance \(\frac{\nu_i}{\nu_i-2}u_i^2\), as mc.margins = "S1" in propagate. For a sampler or a sample, the distribution is that of its values, and \(x_i\) and \(u_i\) are their mean and standard deviation (df = Inf).

Sampling. \(\mathbf{Z} \sim N(\mathbf{0}, \mathbf{R}_Z)\) is drawn with the copula correlation matrix \(\mathbf{R}_Z\), and \(X_i = Q_i(\Phi(Z_i))\) with the quantile function \(Q_i\) of the margin (Gaussian copula, Sklar's theorem; Possolo 2010). The rank correlations of the samples are those of the copula, \(\rho_S = \frac{6}{\pi}\arcsin\frac{\rho_Z}{2}\), for any margins. For samplers and samples, the nsim values are sorted and assigned in the rank order of \(Z_i\) (rank reordering), which is the same construction with the empirical quantile function: each input keeps exactly the values drawn from its distribution, and only their order is changed. All 34 generators of rDistr, all random number generators of R and the distributions of measured data can thus be correlated.

Correlation. With cor.type = "copula", \(\mathbf{R}_Z = \mathbf{R}\); the Pearson correlation of non-Gaussian margins is then slightly smaller (e.g. 0.786 instead of 0.80 for two rectangular inputs). With cor.type = "pearson", each element of \(\mathbf{R}_Z\) is found by root finding such that \(\mathrm{cor}(X_i, X_j) = R_{ij}\), with the expectation evaluated by Gauss-Hermite quadrature (NORTA, Cario & Nelson 1997); for two rectangular margins this is \(\rho_Z = 2\sin(\pi R_{ij}/6)\). The GUM uncertainty framework and the Monte Carlo method then use the same input covariance matrix, as required for the validation of the GUF by the Monte Carlo method (JCGM 101:2008, 8). Correlations that cannot be attained with the given margins give an error; for t-margins with \(\nu_i \le 4\) (the sample correlation has no finite variance), the copula value \(\mathbf{R}_Z = \mathbf{R}\) is kept.

Use with propagate. The samples carry the attribute "budget" with the stated estimates, standard uncertainties, degrees of freedom and distributions. propagate(..., type = "sim", cov = "full") then uses these values for the GUM uncertainty framework, together with the correlation matrix of the samples, and the samples themselves for the Monte Carlo simulation. The adaptive Monte Carlo procedure (adaptive = TRUE) draws new samples from the budget without correlation and must not be used with corDist samples; increase nsim instead.

References

Copulas:
Copulas for uncertainty analysis.
Possolo A.
Metrologia (2010), 47: 262-271.
tools:::Rd_expr_doi("10.1088/0026-1394/47/3/017")

Modeling and generating random vectors with arbitrary marginal distributions and correlation matrix.
Cario MC, Nelson BL.
Technical report, Northwestern University (1997).

Probability distributions of input quantities:
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")

See Also

propagate, makeDat, rDistr

Examples

## Two correlated rectangular inputs (r = 0.8), Y = X1 + X2:
## u(y) = sqrt(1 + 1 + 2 * 0.8) = 1.897 for GUF and Monte Carlo
R <- matrix(c(1, 0.8, 0.8, 1), 2)
X <- corDist(mean = c(X1 = 10, X2 = 5), u = c(1, 1), dist = "rect", cor = R, nsim = 1E5)
cor(X)[1, 2]
RES <- propagate(expression(X1 + X2), X, type = "sim", cov = "full", check = FALSE)
c(RES$prop["u.1"], RES$sim["u.MC"])

## Mixed margins: t (10 degrees of freedom), curvilinear trapezoidal
## and a lognormal given by a sampler
R3 <- matrix(c(1, 0.5, -0.3, 0.5, 1, 0, -0.3, 0, 1), 3)
X3 <- corDist(mean = c(a = 1, b = 2, z = NA), u = c(0.1, 0.2, NA),
              dist = list("t", "ctrap", function(n) rlnorm(n, 0.5, 0.3)),
              df = c(10, Inf, Inf), d = c(NA, 0.05, NA), cor = R3, nsim = 1E5)
round(cor(X3), 2)
propagate(expression(a * b + z), X3, type = "sim", cov = "full", check = FALSE)

## Internal samplers of the package (trapezoidal, 4P beta) and a measured sample
meas <- c(50.3, 49.1, 51.2, 50.8, 48.7, 50.1, 49.9, 51.5, 50.4, 49.6)
X4 <- corDist(mean = c(A = NA, B = NA, M = NA), u = NA,
              dist = list(function(n) propagate:::rtrap(n, 0, 1, 2, 3),
                          function(n) propagate:::rbeta2(n, 2, 5, 10, 12), meas),
              cor = matrix(c(1, 0.7, 0.4, 0.7, 1, 0.2, 0.4, 0.2, 1), 3), nsim = 1E5)
round(cor(X4), 2)
attr(X4, "budget")