propagateMV: Uncertainty propagation for several output quantities (JCGM 102:2011)

Description

Propagation of uncertainty to several output quantities that depend on common input quantities, by the GUM uncertainty framework (law of propagation of uncertainty for the covariance matrix of the outputs) and by the Monte Carlo method on a common input sample, as described in GUM Supplement 2 (JCGM 102:2011). The function uses propagate internally, so that all its input types and Monte Carlo input models are available.

Usage

propagateMV(expr, data, type = c("stat", "raw", "sim"), cov = c("diag", "full"),
            nsim = 1000000, alpha = 0.05, mc.margins = c("S1", "mvt"), ...)
# S3 method for propagateMV
print(x, digits = 7, ...)
# S3 method for propagateMV
summary(object, ...)
# S3 method for propagateMV
plot(x, pairs = NULL, npoints = 3000, ...)

Value

An object of class propagateMV, a list with

guf: results of the GUM uncertainty framework: estimates y, covariance matrix U, standard uncertainties u, correlation matrix r, coverage factors kp and kq, sensitivity matrix Cx.
mc: Monte Carlo results: y, U, u, r, kp, kq and the coverage intervals ci of the single outputs.
resSIM: the matrix of Monte Carlo output values (trials $\times$ outputs).
datSIM: the joint Monte Carlo input sample.
covMat: the covariance matrix $\mathbf{U}_x$ of the inputs.
expr, names, alpha, nsim, type, mc.margins, df.MC: the model, output names and settings.
runs: the propagate objects of the single outputs (univariate GUF results, sensitivities, expanded uncertainties).

Arguments

expr: an expression with one element per output quantity, preferably named, e.g. expression(R = (V/I) * cos(phi), X = (V/I) * sin(phi)), or a list of expressions. Unnamed outputs are called Y1, Y2, ...
data: the input data as for propagate, containing all variables of all output models.
type: "stat", "raw" or "sim", as in propagate.
cov: "diag", "full" or an external covariance matrix of the inputs, as in propagate.
nsim: the number of Monte Carlo trials.
alpha: 1 - coverage probability of the coverage regions and intervals.
mc.margins: the Monte Carlo input model for Type A inputs, as in propagate.
x, object: an object returned by propagateMV.
digits: the number of significant digits printed.
pairs: a list of pairs of output indices to plot; default all pairs.
npoints: the number of Monte Carlo points plotted per panel.
...: further arguments to propagate (e.g. df.tot, ws.group) or to plot.

Author

Andrej-Nikolai Spiess

Details

Joint Monte Carlo sample. propagate is called once with a carrier model: the first output model plus $0 \cdot x$ for every input $x$ that it does not contain. This generates the joint sample of all input quantities with the chosen input model (Gaussian copula with individual degrees of freedom, multivariate t-distribution, user samples, ...). All output models are evaluated on this same sample, so that the Monte Carlo results describe the joint distribution of the outputs (JCGM 102:2011, 7): estimates $\tilde{\mathbf{y}}$ (means), covariance matrix $\mathbf{U}_{\tilde{y}}$, correlation matrix and, for each output, the probabilistically symmetric coverage interval (JCGM 101:2008, 7.7.1). Trials with non-finite outputs are removed with a message.

GUM uncertainty framework. propagate is called once per output to obtain its sensitivity coefficients at the estimates (symbolic derivatives, or numerical derivatives if necessary). With the sensitivity matrix $\mathbf{C}_x$ (outputs $\times$ inputs) and the covariance matrix $\mathbf{U}_x$ of the inputs, $$\mathbf{U}_y = \mathbf{C}_x \mathbf{U}_x \mathbf{C}_x^\top$$ (JCGM 102:2011, 6.2.1.3, Eq. 3).

Coverage regions for coverage probability $p = 1 - \alpha$ (JCGM 102:2011, 6.5 and 7.7):

hyper-ellipsoidal: $(\boldsymbol{\eta} - \mathbf{y})^\top \mathbf{U}_y^{-1} (\boldsymbol{\eta} - \mathbf{y}) = k_p^2$. GUF: $k_p = \sqrt{\chi^2_{p}(m)}$ for $m$ outputs (Gaussian); Monte Carlo: $k_p^2$ is the $p$-quantile of the squared Mahalanobis distances of the trials from the Monte Carlo mean.
hyper-rectangular: $\mathbf{y} \pm k_q \mathbf{u}(\mathbf{y})$. GUF: $k_q = \Phi^{-1}((1 + p^{1/m})/2)$ (independent Gaussian components, e.g. 2.24 for $m = 2$, $p = 0.95$, JCGM 102:2011, Table 2); Monte Carlo: $k_q$ is the $p$-quantile of $\max_j |y_j - \tilde{y}_j|/u(\tilde{y}_j)$.

$k_p$ is not available if the covariance matrix of the outputs is singular (linearly dependent outputs).

Plot. For each pair of outputs, the Monte Carlo sample with the bivariate $100(1-\alpha)$% ellipses of the GUF (dashed) and the Monte Carlo method (solid), each with $k = \sqrt{\chi^2_{1-\alpha}(2)}$ from the respective estimates and covariance matrix of the pair. The Monte Carlo ellipse therefore uses the Gaussian factor, not the Monte Carlo $k_p$, and is not the Monte Carlo coverage region.

The function reproduces the examples of JCGM 102:2011, clause 9; see demo("JCGM102", package = "propagate").

References

Evaluation of measurement data - Supplement 2 to the "Guide to the expression of uncertainty in measurement" - Extension to any number of output quantities.
JCGM 102:2011. Joint Committee for Guides in Metrology (2011).
tools:::Rd_expr_doi("10.59161/JCGM102-2011")

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")

Examples

Run this code

## JCGM 102:2011, 9.4: resistance, reactance and impedance from six sets of
## simultaneous observations (GUM H.2 plus one set), multivariate t-distribution
## of the inputs with nu = n - N = 3 degrees of freedom.
data(S2.9.4)
EXPR <- expression(R = (V/I) * cos(phi), X = (V/I) * sin(phi), Z = V/I)
RES <- propagateMV(EXPR, S2.9.4, type = "raw", cov = "full", mc.margins = "mvt", nsim = 1E5)
RES
## JCGM 102:2011, Table 11: GUF u = 0.058, 0.241, 0.193 Ohm;
## MCM u = 0.130, 0.536, 0.429 Ohm, r(R, X) = -0.587 (M = 1e6)
plot(RES)

## JCGM 102:2011, 9.2, case 3: common rectangular effect X3,
## input samples from an uncertainty budget
BUD <- data.frame(name = c("X1", "X2", "X3"), value = 0, u = c(1, 1, NA),
                  a = c(NA, NA, 3 * sqrt(3)), dist = c("normal", "normal", "rect"))
EX2 <- expression(Y1 = X1 + X3, Y2 = X2 + X3)
DAT <- makeDat(EX2, data = BUD, output = "sim", nsim = 1E5)
RES2 <- propagateMV(EX2, DAT, type = "sim")
c(RES2$guf$kp, RES2$guf$kq)   ## GUF: 2.45, 2.24
c(RES2$mc$kp, RES2$mc$kq)     ## MCM: 2.28, 1.87 (JCGM 102:2011, Table 5)

Run the code above in your browser using DataLab