Error propagation in R: addition in quadrature and the delta 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.

Not known to all physicists and engineers, however, is that the addition in quadrature procedure is exactly the same as what statisticians call the delta method. But the opposite is also true: not all statisticians seem to realize that the delta method is identical to the addition in quadrature method used by physicists and engineers.

The following R code demonstrates how to calculate the uncertainty in physical measurements using the delta method. The resulting uncertainty is always identical to that obtained by the addition in quadrature procedure.

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 derivation of the delta method (and demonstrations of its application) can be found in Meeker and Escobar’s 1998 book Statistical methods for reliability data.

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)

##example from Taylor, pp. 69-70

#problem: calculate the the refractive index (and its uncertainty) using Snell's law

#abbreviations used:
#i = angle of incidence
#r = angle of refraction
#n = refractive index

#the refractive index is calculated as n = sin(i)/sin(r)

#following results were obtained when measuring angles i and r:
i <- 20 #mean of i (in degrees)
sdi <- 1 #standard deviation of i
r <- 13 #mean of r (in degrees)
sdr <- 1 #standard deviation of r

#convert degrees to radians (by multiplying with pi/180)
i <- i*(pi/180); sdi <- sdi*(pi/180); r <- r*(pi/180); sdr <- sdr*(pi/180)

#construct variance covariance matrix of i and r
vcovData <- matrix(nrow=2, ncol=2)
vcovData[1,1] <- sdi^2 #variance = (standard deviation)^2
vcovData[2,2] <- sdr^2
#measures of i and r were assumed to be independent (i.e., cov[i,r]=0)
vcovData[1,2] <- vcovData[2,1] <- 0
colnames(vcovData) <- rownames(vcovData) <- c("i", "r")
vcovData

#calculate uncertainty with delta method

#partial derivative of sin(i)/sin(r) with respect to i
dgdi <- deriv(~ sin(i)/sin(r), "i")
#partial derivative of sin(i)/sin(r) with respect to r
dgdr <- deriv(~ sin(i)/sin(r), "r")
#evaluate these partial derivatives at the means of i and r, respectively
dgdir <- c(attr(eval(dgdi),"gradient"), attr(eval(dgdr),"gradient"))
#calculate standard error of refractive index (n)
SEn <- sqrt(t(dgdir)%*%vcovData%*%dgdir)

#results:
#estimate of refractive index (report 3 significant figures)
(n <- round(sin(i)/sin(r), 2))
#estimate of uncertainty (or error) in refractive index
round(SEn, 2)
#fractional uncertainty
round(SEn/n, 2)
#these results are identical to those reported by Taylor who used the
#addition in quadrature procedure


#applying the deltaMethod function from the car package

#construct variance covariance matrix for i and r
vcovData2 <- vcovData
colnames(vcovData2) <- rownames(vcovData2) <- c("angleI", "angleR")
vcovData2

#construct vector containing the means of i and r
meansMeasures <- c(i, r)
names(meansMeasures) <- c("angleI", "angleR")

#calculate uncertainty with delta method
err <- deltaMethod(meansMeasures, g="sin(angleI)/sin(angleR)", vcov.=vcovData2)
err
#the R package "propagate" also implements the delta method
#see the function called propagate in that R package



##example from Taylor, pp. 68-69

#problem: calculate the gravitational acceleration (and its uncertainty)
#using a simple pendulum

#abbreviations used:
#g = gravitational acceleration
#l = length of pendulum
#T = period of pendulum

#note that the period of a pendulum is given by T=2*pi*sqrt(l/g)
#thus, g=4*(pi^2)*(l/T^2)

#following results were obtained when measuring l and T:
meanLength <- 92.95 #mean l (in cm)
sdLength <- .1 #standard deviation of l
meanTime <- 1.936 #mean T (in seconds)
sdTime <- .004 #standard deviation of T

#construct variance covariance matrix of l and T
vcovData <- matrix(nrow=2, ncol=2)
vcovData[1,1] <- sdLength^2 #variance = (standard deviation)^2
vcovData[2,2] <- sdTime^2
#measures of l and T were assumed to be independent (i.e., cov[l,T]=0)
vcovData[1,2] <- vcovData[2,1] <- 0
colnames(vcovData) <- rownames(vcovData) <- c("length", "time")
vcovData

#construct vector containing the means of length and period
meansMeasures <- c(meanLength, meanTime)
names(meansMeasures) <- c("length", "time")
meansMeasures

#calculate uncertainty using the delta method
err <- deltaMethod(meansMeasures, g="4*(pi)^2*length/(time^2)", vcov.=vcovData)
err

#estimate of g (report 3 significant figures)
round(err$Estimate)
#estimate of uncertainty (or error) in g
round(err$SE)
#note: estimates are identical to those reported by Taylor who used the
#addition in quadrature procedure



##example from Taylor, pp. 213-214

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

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

meanAlpha <- mean(alpha) #mean alpha
#note: Taylor used the population standard deviation in this example
#the population standard deviation can be obtained in R as follows:
#var(x)*(n-1)/n, where n is the number of observations of variable x
varAlpha <- var(alpha)*(length(alpha) - 1)/length(alpha) #variance alpha
meanBeta <- mean(beta) #mean beta
varBeta <- var(beta)*(length(beta) - 1)/length(beta) #variance beta
covAlphaBeta <- cov(alpha,beta)*(length(alpha) - 1)/length(alpha) #cov(alpha,beta)

#construct variance covariance matrix of alpha and beta
vcovData <- matrix(nrow=2, ncol=2)
vcovData[1,1] <- varAlpha
vcovData[2,2] <- varBeta
vcovData[1,2] <- vcovData[2,1] <- covAlphaBeta
colnames(vcovData) <- rownames(vcovData) <- c("alpha", "beta")
vcovData

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

#calculate uncertainty with delta method
err <- deltaMethod(meansMeasures, g="alpha+beta", vcov.=vcovData)
err

#estimate of g (report 2 significant figures)
round(err$Estimate)
#estimate of uncertainty (or error) in alpha+beta
round(err$SE)
#note: estimates are identical to those reported by Taylor who
#used the addition in quadrature procedure