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, "to", r, "\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, "to", r, "\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, "to", r, "\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.