Error propagation in R: addition in quadrature versus the bootstrap method

Physicists and engineers often have to calculate the uncertainty in a derived quantity.

For instance, a test engineer measures two angles. The uncertainty (or error) in these measurements appears to be +/- 1 degree. Subsequently, the engineer calculates the sum of these two angles. This sum is a derived quantity. But note that this derived quantity is composed of two measurements each having their own uncertainties, so what is the uncertainty (or error) in this derived quantity? In other words, how propagates the uncertainty from the measured quantities (the two angles) to a derived quantity (sum of two angles). For calculating the uncertainty in this derived quantity, physicists and engineers rely on what is called the addition in quadrature procedure.

The addition in quadrature procedure provides an estimate of the Standard Deviation Of the Mean (SDOM), which is a quantification of the uncertainty in a derived quantity. This SDOM, in turn, can be used for constructing confidence intervals for the derived quantity.
Bootstrap methods (such as the bias-corrected bootstrap method) provide another way of obtaining confidence intervals for the derived quantity.

The following R code shows how to compute confidence intervals for a derived physical quantity with (1) the addition in quadrature procedure, and (2) the bias-corrected bootstrap method. The code demonstrates that these two methods usually yield very similar confidence intervals.

Also note that the bootstrap confidence intervals have a coverage probability that is usually close (or equal) to the nominal coverage probability. That is, when the nominal coverage probability is 95% (i.e., a 95% confidence interval), the coverage probability of the bootstrap intervals seems to be consistently close (or equal) to this 95%.

Towards the end I added
R code for running simulations that will calculate the coverage probabilities for both the addition in quadrature procedure and the bootstrap method.
I encourage you to experiment with these simulations for calculating the coverage probability. But you need to have some patience with these simulations, since running them may take some time…

A derivation of the addition in quadrature procedure (and lots of demonstrations of its application) can be found in John Taylor’s 1997 book An introduction to error analysis: the study of uncertainties in physical measurements.
A step-by-step explanation of the bias-corrected bootstrap method is given by Manly in his 2007 book Randomization, bootstrap and Monte Carlo methods in biology.

Suggestions and/or questions? Please contact Stefan Gelissen (email: info at datall-analyse.nl).
See this page for an overview of all of Stefan’s R code blog posts.

R code

library(car)
library(mvtnorm)

##independent measurements

#example with independent measurements from Taylor, pp. 104-105

#problem: estimate the area of a rectangular plate and calculate the
#uncertainty in this area estimate

#abbreviations used:
#l = length of plate
#b = breath of plate

#following results were obtained when measuring l and b (in cm)
l <- c(24.25,24.26,24.22,24.28,24.24,24.25,24.22,24.26,24.23,24.24)
b <- c(50.36,50.35,50.41,50.37,50.36,50.32,50.39,50.38,50.36,50.38)

meanL <- mean(l) #mean l
meanB <- mean(b) #mean b
nl <- length(l) #number of l measurements
nb <- length(b) #number of b measurements



##addition in quadrature procedure
#independent measurements

#note: assume that length and breath both follow a normal distribution

#construct variance covariance matrix of l and b
vcovData <- matrix(nrow=2, ncol=2)
vcovData[1,1] <- var(l)/nl
vcovData[2,2] <- var(b)/nb
#measures of l and b are assumed to be independent (i.e., cov[l,b]=0)
vcovData[1,2] <- vcovData[2,1] <- 0
colnames(vcovData) <- rownames(vcovData) <- c("length", "breadth")
round(vcovData,6)

#construct vector containing the means of l and b
meansMeasures <- c(meanL, meanB)
names(meansMeasures) <- c("length", "breath")
meansMeasures

#calculate uncertainty with delta method
#note: the delta method and addition in quadrature procedure are identical
#for more information on this identity, see: http://blogs2.datall-analyse.nl/2016/02/18/error_propagation_r_addition_in_quadrature/
err <- deltaMethod(meansMeasures, g="length*breath", vcov.=vcovData)
err

#obtain 95% confidence interval
alpha <- .05
#determine effective degrees of freedom with Welch-Satterthwaite formula
ul <- sd(l)/sqrt(nl) #l's standard deviation of the mean
l.df <- nl-1 #degrees of freedom associated with l
df.dl <- deriv(~meanL*meanB, "meanL") #partial derivative of l*b with respect to l
cl <- attr(eval(df.dl),"gradient") #evaluate the partial derivative
ub <- sd(b)/sqrt(nb) #b's standard deviation of the mean
b.df <- nb-1 #degrees of freedom associated with b
df.db <- deriv(~meanL*meanB, "meanB") #partial derivative of l*b with respect to b
cb <- attr(eval(df.db),"gradient") #evaluate the partial derivative
#Welch-Satterthwaite
ws.df <- floor(err$SE^4/sum(c( (cl*ul)^4/l.df, (cb*ub)^4/b.df)))
#confidence interval
err$Estimate + c(1,-1)*qt(alpha/2, df=ws.df)*err$SE



##bias-corrected bootstrap method
#independent measurements

#equation to be evaluated
equation <- function(x, y) x*y
hatXY <- equation(meanL, meanB)

#initialize bootstrap
B <- 4000 #number of bootstrap samples
alpha <- .05 #confidence level
estimateB <- numeric(B)

meanX <- meanL
sdX <- sd(l)
meanY <- meanB
sdY <- sd(b)
nx <- nl
ny <- nb
#number of paired (or simultaneous) observations of l and b
n <- nrow(cbind(l,b))

#perform bootstrap
for (i in 1:B) {
  xb <- meanX + rt(n, df=nx-1)*sdX
  yb <- meanY + rt(n, df=ny-1)*sdY
  mxb <- mean(xb)
  myb <- mean(yb)
  estimateB[i] <- equation(x=mxb, y=myb)
}

#check results bootstrap
hist(estimateB)
mean(estimateB)
sd(estimateB)

#bias-corrected percentile confidence interval (two sided)
q0 <- mean(estimateB <= hatXY)
zl.bias2 <- round(B*pnorm(2*qnorm(q0)+qnorm(alpha/2)))
zu.bias2 <- round(B*pnorm(2*qnorm(q0)+qnorm(1-alpha/2)))

c(sort(estimateB)[zl.bias2], sort(estimateB)[zu.bias2])

#standard (or simple) percentile confidence intervals
quantile(estimateB, probs=c(alpha/2, 1-alpha/2))



##Monte Carlo simulation
#compare bootstrap confidence interval with that obtained by a Monte Carlo simulation
lSim <- mean(l) + rt(10000, df=nl-1)*sd(l)/sqrt(nl)
bSim <- mean(b) + rt(10000, df=nb-1)*sd(b)/sqrt(nb)
areaSim <- lSim*bSim
#results Monte Carlo simulation
hist(areaSim)
mean(areaSim)
sd(areaSim)
#95% confidence interval
quantile(areaSim, probs=c(alpha/2, 1-alpha/2))





##dependent measurements

#example with dependent (correlated) measurements from Taylor, pp. 213-214

#problem: add two angles and calculate the uncertainty in this sum

#following results were obtained when measuring two angles beta and gamma
beta <- c(35,31,33,32,34)
gamma <- c(50,55,51,53,51)

meanBeta <- mean(beta)
varBeta <- var(beta)
meanGamma <- mean(gamma)
varGamma <- var(gamma)
covBetaGamma <- cov(beta,gamma)

nb <- length(beta) #number of beta measurements
ng <- length(gamma) #number of gamma measurements
#number of paired (or simultaneous) observations of beta and gamma
nbg <- nrow(cbind(beta,gamma))



##addition in quadrature procedure
#dependent measurements

#note: assume that beta and gamma follow a multivariate normal distribution

#construct variance covariance matrix of alpha and beta
vcovData <- matrix(nrow=2, ncol=2)
vcovData[1,1] <- varBeta/nb
vcovData[2,2] <- varGamma/ng
vcovData[1,2] <- vcovData[2,1] <- covBetaGamma/nbg
colnames(vcovData) <- rownames(vcovData) <- c("beta", "gamma")
vcovData

#construct vector containing the means of alpha and beta
meansMeasures <- c(meanBeta, meanGamma)
names(meansMeasures) <- c("beta", "gamma")
meansMeasures

#calculate uncertainty with delta method (=addition in quadrature procedure)
err <- deltaMethod(meansMeasures, g="beta+gamma", vcov.=vcovData)
err

#95% confidence intervals are not available for the addition in quadrature procedure:
#note that the Welch-Satterthwaite formula for calculating the degrees of
#freedom is only valid when measurements are independent, which was not the case
#with the measurements of beta and gamma
#however, the bias-corrected bootstrap method yields valid confidence intervals
#in case of dependent measurements
#simulation results (reported at the end of this R code) demonstrate that
#these bootstrap confidence intervals have a coverage probability that is close
#to the nominal coverage probability



##bias-corrected bootstrap method
#dependent measurements

#equation to be evaluated
equation <- function(x, y) x+y
hatXY <- equation(meanBeta, meanGamma)

#initialize bootstrap
B <- 4000 #number of bootstrap samples
alpha <- .05
estimateB <- numeric(B)

n <- nrow(cbind(beta,gamma))
meanX <- meanBeta
varX <- varBeta
meanY <- meanGamma
varY <- varGamma
#construct variance covariance matrix
vcovXY <- matrix(nrow=2, ncol=2)
vcovXY[1,1] <- varBeta
vcovXY[2,2] <- varGamma
vcovXY[1,2] <- vcovXY[2,1] <- covBetaGamma

#perform bootstrap
for (i in 1:B) {
  mus <- rep(c(meanX, meanY), each=n)
  #for the multivariate t-distribution, we will assume that each
  #of the marginal t-distributions has n-1 degrees of freedom
  xyb <- mus + rmvt(n, sigma=vcovXY, df=n-1, method="chol")
  mxb <- mean(xyb[,1])
  myb <- mean(xyb[,2])
  estimateB[i] <- equation(x=mxb, y=myb)
}

#check results bootstrap
hist(estimateB)
mean(estimateB)
sd(estimateB)

#bias-corrected percentile confidence interval (two sided)
q0 <- mean(estimateB <= hatXY)
zl.bias2 <- round(B*pnorm(2*qnorm(q0)+qnorm(alpha/2)))
zu.bias2 <- round(B*pnorm(2*qnorm(q0)+qnorm(1-alpha/2)))

c(sort(estimateB)[zl.bias2], sort(estimateB)[zu.bias2])

#standard (or simple) percentile confidence intervals
quantile(estimateB, probs=c(alpha/2, 1-alpha/2))



##Monte Carlo simulation
#compare bootstrap confidence interval with that obtained by a Monte Carlo simulation
mus <- rep(c(mean(beta), mean(gamma)), each=10000)
bgSim <- mus + rmvt(10000, sigma=vcovData, df=nbg-1, method="chol")
sumSim <- bgSim[,1]+bgSim[,2]
#results Monte Carlo simulation
hist(sumSim)
mean(sumSim)
sd(sumSim)
#95% confidence interval
quantile(sumSim, probs=c(alpha/2, 1-alpha/2))



##Monte Carlo simulation using a Gaussian copula with Student margins
library(copula)

#obtain start values for fitting the Gaussian copula
sdBeta <- sqrt((nb-1)/nb)*sd(beta)
sdGamma <- sqrt((ng-1)/ng)*sd(gamma)
corBetaGamma <- sin(cor(beta, gamma, method="kendall")*pi/2)
start <- c(meanBeta, sdBeta, meanGamma, sdGamma, corBetaGamma)

#fit Gaussian copula
myGaussian <- mvdc(copula=normalCopula(param=1, dispstr="un"),
                   margins=c("norm", "norm"),
                   paramMargins=list(list(mean=0, sd=1), list(mean=0, sd=1)))
fitGaussian <- fitMvdc(data=cbind(beta, gamma), mvdc=myGaussian, start=start)
fitGaussian

#note that the MLE for the standard deviations are biased
summary(lm(beta~1))$sigma #unbiased OLS estimate for the sd of beta
fitGaussian@estimate[2] #biased MLE for the sd of beta
sqrt(nb/(nb-1))*fitGaussian@estimate[2] #unbiased MLE for the sd of beta
summary(lm(gamma~1))$sigma #unbiased OLS estimate for the sd of gamma
fitGaussian@estimate[4] #biased MLE for the sd of gamma
sqrt(ng/(ng-1))*fitGaussian@estimate[4] #unbiased MLE for the sd of gamma

#set up Gaussian copula with Student margins for Monte Carlo simulation
corBetaGammaFit <- fitGaussian@mvdc@copula@parameters
meanBfit <- fitGaussian@estimate[1]
sdBfit <- sqrt(nb/(nb-1))*fitGaussian@estimate[2]
meanGfit <- fitGaussian@estimate[3]
sdGfit <- sqrt(ng/(ng-1))*fitGaussian@estimate[4]

myStudent <- mvdc(copula=normalCopula(param=corBetaGammaFit, dispstr="un", dim=2),
                  margins=c("t", "t"), paramMargins=list(list(df=nb-1), list(df=ng-1)))

#Monte Carlo simulation
bgSim <- rMvdc(10000, myStudent)
bgSim[,1] <- meanBfit + bgSim[,1]*sdBfit/sqrt(nb)
bgSim[,2] <- meanGfit + bgSim[,2]*sdGfit/sqrt(ng)
sumSim <- bgSim[,1]+bgSim[,2]
#results Monte Carlo simulation
hist(sumSim)
mean(sumSim)
sd(sumSim)
#95% confidence interval
quantile(sumSim, probs=c(alpha/2, 1-alpha/2))

#shortest 95% coverage interval
nsim <- length(sumSim)
sortedSim <- sort(sumSim)
covInt <- sapply(1:(nsim-round((1-alpha)*nsim)), function(i) {
  sortedSim[1+round((1-alpha)*nsim)+(i-1)]-sortedSim[1+(i-1)]})
#lower limit
(lcl <- sortedSim[which(covInt==min(covInt))])
#upper limit
(ucl <- sortedSim[1+round((1-alpha)*nsim)+(which(covInt==min(covInt))-1)])
#draw limits
hist(sumSim)
abline(v=lcl, col="blue", lty=2)
abline(v=ucl, col="blue", lty=2)





##estimate coverage probability of confidence interval
#independent measurements

#true means and sd's (mimicking rectangular plate data from Taylor, pp. 104-105)
meanLT <- 24; sdLT <- .02 #assume that length follows a normal distribution
meanBT <- 50; sdBT <- .02 #assume that breath follows a normal distribution
XYt <- meanLT*meanBT

#initialize procedure for estimating coverage probability
nsim <- 10000 #number of simulations
alpha <- .05 #nominal coverage probability is equal to 1-alpha
m <- 20 #sample size

##coverage probability of delta method (=addition in quadrature)
vcovB <- matrix(nrow=2, ncol=2)
uclB <- numeric(nsim)
lclB <- numeric(nsim)

for (j in 1:nsim) {
  xb <- rnorm(m, meanLT, sdLT)
  yb <- rnorm(m, meanBT, sdBT)
  mxb <- mean(xb)
  myb <- mean(yb)
  vcovB[1,1] <- var(xb)/m
  vcovB[2,2] <- var(yb)/m
  vcovB[1,2] <- vcovB[2,1] <- 0
  colnames(vcovB) <- rownames(vcovB) <- c("xb", "yb")
  mVec <- c(mxb,myb)
  names(mVec) <- c("xb", "yb")
  errB <- deltaMethod(mVec, g="xb*yb", vcov.=vcovB)
  sdom <- errB$SE
  ux <- sd(xb)/sqrt(m)
  uy <- sd(yb)/sqrt(m)
  ws.df <- floor(sdom^4/sum(c( (myb*ux)^4/(m-1), (mxb*uy)^4/(m-1) ) ))
  lclB[j] <- errB$Estimate + qt(alpha/2, df=ws.df)*sdom
  uclB[j] <- errB$Estimate - qt(alpha/2, df=ws.df)*sdom
}

#coverage probability
mean(lclB < XYt & uclB > XYt) #result will be 95% (approximately)
#average width of confidence intervals
mean(uclB-lclB) #result will be 1.01 (approximately)



##coverage probability of bias-corrected bootstrap method
B <- 4000
estimateBB <- numeric(B)
uclB <- numeric(nsim)
lclB <- numeric(nsim)

for (j in 1:nsim) {
  xb <- rnorm(m, meanLT, sdLT)
  yb <- rnorm(m, meanBT, sdBT)
  mxb <- mean(xb)
  myb <- mean(yb)
  sdxb <- sd(xb)
  sdyb <- sd(yb)
  hatXYb <- mxb*myb
  
  for (i in 1:B) {
    xbB <- mxb + rt(m, df=m-1)*sdxb
    ybB <- myb + rt(m, df=m-1)*sdyb
    mxbB <- mean(xbB)
    mybB <- mean(ybB)
    estimateBB[i] <- mxbB*mybB}
  
  q0 <- mean(estimateBB <= hatXYb)
  lclB[j] <- sort(estimateBB)[round(B*pnorm(2*qnorm(q0)+qnorm(alpha/2)))]
  uclB[j] <- sort(estimateBB)[round(B*pnorm(2*qnorm(q0)+qnorm(1-alpha/2)))]
}

#coverage probability
mean(lclB < XYt & uclB > XYt) #result will be 95% (approximately)
#average width of confidence intervals
mean(uclB-lclB) #result will be 1.02 (approximately)



##coverage probability of Monte Carlo simulation
nMC <- 10000
uclB <- numeric(nsim)
lclB <- numeric(nsim)

for (j in 1:nsim) {
  xb <- rnorm(m, meanLT, sdLT)
  yb <- rnorm(m, meanBT, sdBT)
  mxb <- mean(xb)
  myb <- mean(yb)
  sdxb <- sd(xb)
  sdyb <- sd(yb)
  
  xSim <- mxb + rt(nMC, df=m-1)*sdxb/sqrt(m)
  ySim <- myb + rt(nMC, df=m-1)*sdyb/sqrt(m)
  hatXYb <- xSim*ySim
  
  lclB[j] <- quantile(hatXYb, probs=c(alpha/2))
  uclB[j] <- quantile(hatXYb, probs=c(1-alpha/2))
}

#coverage probability
mean(lclB < XYt & uclB > XYt) #result will be 95% (approximately)
#average width of confidence intervals
mean(uclB-lclB) #result will be 1.03 (approximately)



##summary simulation results##

#1) addition in quadrature procedure:
#coverage probability=95%; average width of intervals=1.01

#2) bootstrap method:
#coverage probability=95%; average width of intervals=1.02

#3) Monte Carlo simulation:
#coverage probability=95%; average width of intervals=1.03


##results of identical simulation, but with sample size m=10 (instead of m=20)

#1) addition in quadrature procedure:
#coverage probability=95%; average width of intervals=1.50

#2) bootstrap method:
#coverage probability=96%; average width of intervals=1.53

#3) Monte Carlo simulation:
#coverage probability=96%; average width of intervals=1.55


##results of identical simulation, but with sample size m=5

#1) addition in quadrature procedure:
#coverage probability=95%; average width of intervals=2.45

#2) bootstrap method:
#coverage probability=97%; average width of intervals=2.65

#3) Monte Carlo simulation:
#coverage probability=97%; average width of intervals=2.64


##results of simulations with dependent measures

#initialization of simulations:

#true means and sd's (mimicking two angles data from Taylor, pp. 213-214)
#assume that the two angles follow a multivariate normal distribution, specified by:
#meanBetaT=30; sdBetaT=1.5; meanGammaT=50; sdGammaT=2; covBetaGammaT=-3
#derived quantity: XYt=meanBetaT+meanGammaT

#number of simulations / nominal coverage probability / sample size
#nsim=10000; alpha=.05; m=5

#1) bootstrap method:
#coverage probability=95%; average width of intervals=1.18

#2) Monte Carlo simulation:
#coverage probability=95%; average width of intervals=1.17