uncertMC {metRology}

Monte Carlo evaluation of measurement uncertainty.

Description

uncertMC estimates measurement uncertainty from a function, expression or formula by Monte Carlo simulation.

Usage

uncertMC(expr, x, u, method = "MC", df, cor, cov, distrib, distrib.pars, 
   B = 200, keep.x = TRUE, vectorized=TRUE, ...)

Arguments

expr

An expression, function, or formula with no left-hand side (e.g. ~a*x+b*x^2) which can be evaluated in the environment x to provide a numeric value.

x

A named list or vector of parameters supplied to expr.

u

A named list or named vector of length length(x) of standard uncertainties.

method

Method of uncertainty evaluation. The only method currently supported by uncertMC is "MC". If any other method is specified, control is passed to uncert.

df

A named list or named vector of degrees of freedom. df can be a partial named list if not all distributions (see below) use degrees of freedom.

cor, cov

Optional (square, symmetric) correlation or covariance matrices, respectively. If neither is specified, uncertMC assumes independent variables.

distrib

A character vector of length length(x) or a named list of names of distribution functions associated with u. See Details for defaults.

distrib.pars

A named list of lists of parameters describing the distributions associated with u to be passed to the relevant distribution function. If distrib is present but distrib.pars is not, an attempt is made to set defaults based on other parameters; see Details.

B

Number of Monte Carlo replicates.

keep.x

If TRUE, the simulated replicates of x are included in the return object.

vectorized

If TRUE, expr is assumed to take vector arguments. If FALSE, expr is treated as if it takes scalar arguments. See Details for the difference.

...

Additional parameters to be passed to a function (for the function method) or used in an expression (for expression or formula method).

Details

Although most likely to be called by uncert, uncertMC may be called directly.

If any of x, u, df, distrib or distrib.pars are not lists, they are coerced to lists. If x is not named, arbitrary names of the form 'Xn' are applied. If u, df, distrib or distrib.pars do not have names, the names will be set to names(x) if they are of length exactly length(x); if not, an error is returned.

For Monte Carlo evaluation, distributions and distribution parameters are needed but defaults are used if some or all are absent. If distrib is missing, or if it is a list with some members missing, the distribution is assumed Normal and any missing member of distrib is set to "norm".

Distributions are usually identified by the root of the distribution function name; for example to specify the Normal, distrib$name="norm". At present, only the random value generator (e.g. rnorm) is used. Names of user-specified distributions functions can also be used, provided they have a random value generator named r<dist> where <dist> is the abbreviated distribution. Parameters are passed to distribution functions using do.call, so the function must accept the parameters supplied in distrib.pars.

If distrib.pars or members of it are missing, an attempt is made to deduce appropriate distribution parameters from x, u, df and distrib. In doing so, the following assumptions and values apply for the respective distributions:

norm

mean=x$name, sd=u$name.

unif

min=x-sqrt(3)*u, max=x+sqrt(3)*u.

tri

min=x-sqrt(6)*u, max=x+sqrt(6)*u, mode=x.

t, t.scaled

df=df, mean=x, sd=u.

If either cor or cov are present, a test is made to see if off-diagonal elements are significant. If not, uncertMC treats the values as independent. The test simply checks whether the sum of off-diagonal elements of cor (calculated from cov if cov is present) is bigger than .Machine.double.eps*nrow^2.

Correlation is supported as long as all correlated variables are normally distributed. If correlation is present, uncertMC follows a two-stage simulation procedure. First, variables showing correlation are identified. Following a check that their associated distrib values are all "norm", mvrnorm from the MASS library is called to generate the simulated x values for those variables. Second, any remaining (i.e. independent) variables are simulated from their respective distrib and distrib.pars.

Vectorisation makes a difference to execution speed. If vectorize=TRUE, MC evaluation uses eval using the simulated data as the evaluation environment; if not, apply is used row-wise on the simulated input matrix. This makes an appreciable difference to execution speed (typically eval is faster by a factor of 5 or more) so the default assumes vectorised expressions. However, not all functions and expressions take vector arguments, especially user functions involving complicated arithmetic or numerical solutions. Use vectorize=FALSE for functions or expressions that do not take vector arguments. Note: One common symptom of an expression that does not take vector arguments is an R warning indicating that only the first element (typically of a parameter in x) is used. uncertMC may also return NA for u on attempting to take the sd of a single simulated point.

Value

An object of class uncertMC. See uncertMC-class for details.

Author(s)

S. L. R. Ellison [e-mail]

References

JCGM 100 (2008) Evaluation of measurement data - Guide to the expression of uncertainty in measurement. http://www.bipm.org/utils/common/documents/jcgm/JCGM_100_2008_E.pdf. (JCGM 100:2008 is a public domain copy of ISO/IEC Guide to the expression of uncertainty in measurement (1995) ).

Kragten, J. (1994) Calculating standard deviations and confidence intervals with a universally applicable spreadsheet technique, Analyst, 119, 2161-2166.

Ellison, S. L. R. (2005) Including correlation effects in an improved spreadsheet calculation of combined standard uncertainties, Accred. Qual. Assur. 10, 338-343.

See Also

uncert, uncert-class, uncertMC-class

Examples

  expr <- expression(a+b*2+c*3+d/2)
  x <- list(a=1, b=3, c=2, d=11)
  u <- lapply(x, function(x) x/10)
  u.MC<-uncertMC(expr, x, u, distrib=rep("norm", 4), method="MC")
  print(u.MC, simplify=FALSE)

  #An example with correlation
  u.cor<-diag(1,4)
  u.cor[3,4]<-u.cor[4,3]<-0.5
  u.formc.MC<-uncertMC(~a+b*2+c*3+d/2, x, u, cor=u.cor, keep.x=TRUE)
  u.formc.MC

  #A non-linear example
  expr <- expression(a/(b-c))
  x <- list(a=1, b=3, c=2)
  u <- lapply(x, function(x) x/20)
  set.seed(403)
  u.invexpr<-uncertMC(expr, x, u, distrib=rep("norm", 3), B=999, keep.x=TRUE )
  u.invexpr

  #Look at effect of vectorize
  system.time(uncertMC(expr, x, u, distrib=rep("norm", 3), B=9999, keep.x=TRUE ))
  system.time(uncertMC(expr, x, u, distrib=rep("norm", 3), B=9999, keep.x=TRUE, vectorize=FALSE))