A few weeks ago I mentioned that the propagation of errors is a bit tricky. Say we want to predict the acceleration in an Atwood machine. The machine consists of a string extended over a pulley with masses at either end, of masses M and m, with M > m. The acceleration is given by
where g is the acceleration due to gravity, which we assume is known exactly. Let’s set g = 1, so we’ll have
We previously found by analytic methods that if and
, then
. But it’s instructive to do a simulation.
Specifically, fix some large . For
, let
be normally distributed with mean 100 and standard deviation 1; let
be normally distributed with mean 50 and standard deviation 1; and let
. Then the mean and standard deviation of the
are estimates of the expected acceleration and its error.
This is very easy in R:
n = 10^4;
M = rnorm(n, 100, 1);
m = rnorm(n, 50, 1);
a = (M-m)/(M+m);
When I ran this code I got mean(a) = 0.3333875, sd(a) = 0.00993982. Furthermore, the computed values of a are roughly normally distributed, as shown by this histogram and Q-Q plot. (The line on the Q-Q plot passes through the point (0, mean(a)) and has slope sd(a).)
This works even if the errors are not normally distributed. For example, we can draw the simulated data from a uniform distribution with the given mean and standard deviation:
Mu = runif(n, 100-sqrt(3), 100+sqrt(3))
mu = runif(n, 50-sqrt(3), 50+sqrt(3))
au = (Mu-mu)/(Mu+mu)
I got mean(au) = 0.3332557 and sd(au) = 0.009936519. The distribution of the simulated results is a bit unusual-looking:
There’s also a way to compute an approximation to the error of the result using calculus, but simulation is cheap.