uncertainty-package: tools:::Rd_package_title("uncertainty")

Description

tools:::Rd_package_description("uncertainty")

Arguments

Author

tools:::Rd_package_author("uncertainty")

Maintainer: tools:::Rd_package_maintainer("uncertainty")

Details

The DESCRIPTION file: tools:::Rd_package_DESCRIPTION("uncertainty") tools:::Rd_package_indices("uncertainty") Define an "uncertainty budget" object, including all the involved input quantities. Then estimate the output quantity object by defining a measurement model, using the "uncertainty budget" and applying an estimation method. Print or plot the output quantity estimates or create a "summary uncertainty" object to print or plot the uncertainty contributions to the output quantity described in the measurement model.

References

JCGM 200:2012. International vocabulary of metrology—Basic and general concepts and associated terms (VIM)

JCGM 100:2008. Guide to the expression of uncertainty of measurement

JCGM 101:2008. Supplement 1 Propagation of distributions usign a Monte Carlo method

S.L.R. Ellison and A. Williams (2012) EURACHEM/CITAC Guide CG 4. Quantifying Uncertainty in Analytical Measurement, ISBN 978-0-948926-30-3.

Becker, R.A., Chambers, J.M. and Wilks, A.R. (1988) The New S Language. Wadsworth & Brooks/Cole.

Examples

Run this code

require(mvtnorm)

cor.mat <- matrix(c(1, -0.7, -0.7, 1), 2, 2)

u.budget <- uncertaintyBudget(x=list(name = c("x0", "x1"), 
	mean = c(10, 20), u=c(1, 5), unit = c("kg", "kg"), 
	dof = c(10, 10), type = c("A", "A"),
	description = c("measurand mass", "sample mass"),
	label = c("x[0]", "x[1]"), distribution = c("normal", "normal")),
  y = cor.mat)
u.budget

## Gaussian first order estimates
GFO.res <- uncertainty(x = u.budget,
y = list(measurand_name = "ratio.GFO", 
         measurand_label = expression(ratio[GFO]), 
         measurand_model = "x0/x1", 
	 measurand_description = "ratio of masses at 20 degrees celsius",
         method = "GFO", alpha = 0.05))

contr.GFO <- summary(GFO.res)

## Monte Carlo estimates
MC.res <- uncertainty(x = u.budget, 
y = list(measurand_name = "ratio.MC", 
         measurand_label = expression(ratio[MC]),
         measurand_model = "x0/x1", 
	 measurand_description = "ratio of masses at 20 degrees celsius",
         method = "MC", alpha = 0.05, B = 1e5))

contr.MC <- summary(MC.res)

## print the estimates
MC.res
GFO.res

## print the uncertainty summary
contr.MC
contr.GFO

## Displaying both estimated distributions
if (FALSE) {
plot(MC.res, col = 4, xlab = MC.res$measurand_model)
plot(GFO.res, lty = 2, col = 2, add = TRUE)
legend(0.7, 2.5, legend = c("Monte Carlo", "Gaussian First Order"), 
lty = c(1, 2), col = c(4, 2), lwd = 2, bg = "white")
}

## Display both uncertainty summaries

if (FALSE) {
barplot(cbind(contr.GFO$budget$contrib, contr.MC$budget$contrib), 
beside = TRUE, horiz = TRUE, main = "Uncertainty contribution by method",
xlab = "percent Variance",
names.arg = c(GFO.res$measurand_label, MC.res$measurand_label))
}

##########################
## Example H.1 from GUM ##
##########################

# define the uncertainty budget

u.budget <- uncertaintyBudget(
  x = list(
    name = c("lambda.s", "alpha.s", "theta.bar", "Delta", "delta.alpha",
    "delta.theta", "d.bar", "d.cr",	"d.cnr"),
    label = c("lambda[s]", "alpha[s]", "bar(theta)", "Delta", "delta[alpha]",
    "delta[theta]", "bar(d)", "d[cr]", "d[cnr]"),
    description = c("certified length of the reference gauge block", "reference gauge block cte",
      "mean temperature", "Delta", "difference in the cte", "difference in temperature", 
      "mean difference in length", "d.cr", "d.dnr"),
    mean = c(50.000623,11.5e-6,-1e-1, 0, 0, 0, 2.15e-4, 0, 0),
    unit = c("mm", "oC^-1","oC","oC", "oC^-1", "oC", "mm", "mm", "mm"),
    u = c(25e-6, 1.2e-6, 0.2, 0.35, 0.58e-6, 0.029, 5.8e-6, 3.9e-6, 6.7e-6),
    distribution = c("t", "unif", "unif", "arcsine", "unif", "unif", "t", "t",
      "t"),
    dof = c(18, 1, 1, 1, 50, 2, 24, 5, 8),
    type = c("B", "B", "A", "B", "A", "A", "A", "A", "A")
  ),
  y = diag(1, 9)
)

# define the measurand
measurand_name <- "lambda"
measurand_label <- "lambda"
measurand_model <- paste("(lambda.s*(1+alpha.s*(theta.bar+Delta+delta.theta))",
"+d.bar+d.cr+d.cnr)/(1+(alpha.s+delta.alpha)*(theta.bar+Delta))", sep = "")
measurand_description = "length of end gauge block under calibration at 20 degrees celsius"

# estimate the measurand using the Gaussian First Order method (GUM)

u.GFO <- uncertainty(
	x = u.budget, 
	y = list(measurand_name = measurand_name, 
		measurand_label = measurand_label, 
		measurand_model = measurand_model,
		measurand_description = measurand_description,
		alpha = 0.01,
		method = "GFO"
	)
)

u.GFO
# same result as reported in Table H.1

# estimate the measurand using the Gaussian Second Order method

u.GSO <- uncertainty(
	x = u.budget, 
	y = list(measurand_name = measurand_name, 
		measurand_label = measurand_label, 
		measurand_model = measurand_model,
		measurand_description = measurand_description,
		alpha = 0.01,
		method = "GSO"
	)
) 

u.GSO
# same results as reported in section H.1.6, U(99) = 93 nm,
# the difference is due to rounding error.
# u = 34 nm.

# estimate the measurand using the Monte Carlo method (GUM supplement 1)

u.MC <- uncertainty(
	x = u.budget, 
	y = list(measurand_name = measurand_name, 
		measurand_label = measurand_label, 
		measurand_model = measurand_model,
		measurand_description = measurand_description,
		alpha = 0.01,
		method = "MC", B = 1e5
	)
) 

u.MC
# this result is not reported in the GUM