Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

# Introduction

Error propagation can be a fair bit of work with a calculator, but it’s easy with R. Just use R in repeated calculation with perturbed quantities, and inspect the range of results. Two methods are shown below for inspecting the range: sd() and quantile(), the latter using the fact that area under a normal distribution is 0.68 when calculated between -1 and 1 standard deviation.

# Tests

## Case 1: scale factor

In this case, the answer is simple. If A has uncertainty equal to 0.1, then 10A has uncertainty 1.0.

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 set.seed(123) n <- 500 result <- vector("double", n) A <- 1 Au <- 0.1 # uncertainty in A for (i in 1:n) { Ap <- A + Au * rnorm(n = 1) result[i] = 10 * Ap } D <- 0.5 * (1 - 0.68) r <- quantile(result, probs = c(D, 1 - D)) cat("value:", mean(result), "uncertainty:", sd(result), " range:", r[1], "to", r[2], "\n") 
## value: 10.03 uncertainty: 0.9728  range: 9.047 to 11.02

 1 hist(result) 

The graph indicates that the values are symmetric, which makes sense for a linear operation.

## Case 2: squaring

Here, we expect an uncertainty of 20 percent.

  1 2 3 4 5 6 7 8 9 10 11 12 13 set.seed(123) n <- 500 result <- vector("double", n) A <- 1 Au <- 0.1 # uncertainty in A for (i in 1:n) { Ap <- A + Au * rnorm(n = 1) result[i] = Ap^2 } D <- 0.5 * (1 - 0.68) r <- quantile(result, probs = c(D, 1 - D)) cat("value:", mean(result), "uncertainty:", sd(result), " range:", r[1], "to", r[2], "\n") 
## value: 1.016 uncertainty: 0.1965  range: 0.8184 to 1.213

 1 hist(result) 

## Case 3: a nonlinear function

Here, we have a hyperbolic tangent, so the expected error bar will be trickier analytically, but of course the R method remains simple. (Note that the uncertainty is increased to ensure a nonlinear range of hyperbolic tangent.)

  1 2 3 4 5 6 7 8 9 10 11 12 13 set.seed(123) n <- 500 result <- vector("double", n) A <- 1 Au <- 0.5 # uncertainty in A for (i in 1:n) { Ap <- A + Au * rnorm(n = 1) result[i] = tanh(Ap) } D <- 0.5 * (1 - 0.68) r <- quantile(result, probs = c(D, 1 - D)) cat("value:", mean(result), "uncertainty:", sd(result), " range:", r[1], "to", r[2], "\n") 
## value: 0.7009 uncertainty: 0.233  range: 0.4803 to 0.9065

 1 hist(result) 

# Conclusions

The computation is a simple matter of looping over a perturbed calculation. Here, gaussian errors are assumed, but other distributions could be used (e.g. quantities like kinetic energy that cannot be distributed in a Gaussian manner).

# Further work

1. How large should n be, to get results to some desired resolution?

2. If the function is highly nonlinear, perhaps the mean(result) should be replaced by median(result), or something.