uncertainty-package: Uncertainty Estimation and Contribution Analysis

In uncertainty: Uncertainty Estimation and Contribution Analysis

Uncertainty Estimation and Contribution Analysis

Description

Implements the Gaussian method of first and second order, the Kragten numerical method and the Monte Carlo simulation method for uncertainty estimation and analysis.

Details

The DESCRIPTION file: This package was not yet installed at build time.
Index: This package was not yet installed at build time.
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.

Author(s)

Hugo Gasca-Aragon [aut, cre] (ORCID: <https://orcid.org/0000-0002-5384-2530>)

Maintainer: Hugo Gasca-Aragon <[e-mail]>

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.

See Also

uncertaintyBudget, print.uncertaintyBudget, uncertainty, print.uncertainty, plot.uncertainty, summary.uncertainty, print.summary.uncertainty, plot.summary.uncertainty

Examples

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
## Not run: 
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")

## End(Not run)

## Display both uncertainty summaries

## Not run: 
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))

## End(Not run)

##########################
## 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

uncertainty documentation built on Dec. 15, 2025, 1:06 a.m.