The power and simplicity of propagating errors with Monte Carlo simulations
There's quite a lot about Monte Carlo methods in general at Towards Data Science, but not really much about their very important and useful application to error propagation other than a great introduction by Shuai Guo and a few other articles:
Here, I want to put forward some concrete numerical applications with code for you to actually try and feel first-hand how Monte Carlo methods can be extremely helpful, yet easy-to-implement, for the propagation of errors throughout calculations of almost any kind.
I will begin with a very simple application to propagating the errors during a subtraction operation, to then exemplify how you can use essentially the same idea to propagate errors in virtually any kind of numerical routine from a simple linear regression to a very complex fitting procedure that would be very hard to approach analytically.
Error propagation through Monte Carlo simulations
Error propagation is a fundamental concept in data analysis and scientific computing. When you have measurements with uncertainties, performing mathematical operations on these values will result in propagated errors in the final calculated result. For simple arithmetic operations, error propagation can be done analytically using formulas. If you are interested in analytical error propagation, check out this resource:
However, for more complex operations involving multiple variables and non-linear functions, or for large calculation procedures such as those involved in data fitting or neural network execution, analytical error propagation can quickly become impractical or even impossible.
Monte Carlo simulation offers an alternative approach to error propagation, especially in situations where analytical solutions are challenging or unknown. The Monte Carlo method involves using random sampling to simulate a large number of scenarios, calculating the desired quantity each time, and then analyzing the distribution of results. This statistical approach provides an estimate of the uncertainty in the final result, allowing us to propagate errors through complex calculations.
To demonstrate the application of Monte Carlo simulations for error propagation, let's consider a very simple subtraction operation a – b where both a and b have some associated uncertainties da and db.
(Note: This and all subsequent examples of this article are given in Matlab, but are easily portable to any other scripting language.)
clear; close all; clc
a = 8; da = 2; % Example a = 8 with uncertainty 2
b = 4; db = 1; % Example b = 4 with uncertainty 1
c = a - b; dc = sqrt(da*da + db*db); % Analytical error propagation
disp('Expected')
disp([num2str(c), ' +/- ', num2str(dc)]) % Print
fprintf('n')
% =============== Start here Monte Carlo simulation
c_mc = []; % Will hold all results
for i=1:100000
a_mc = normrnd(a, da); % Get a random number with a Gaussian distribution around a
b_mc = normrnd(b, db); % Same around b
c_mc = [c_mc; a_mc - b_mc]; % Run calculation and store
end
disp('From MonteCarlo with 10 iterations')
disp([num2str(mean(c_mc(1:10))), ' +/- ', num2str(std(c_mc(1:10)))]) % Print result for 10 iterations
fprintf('n')
disp('From MonteCarlo with 100 iterations')
disp([num2str(mean(c_mc(1:100))), ' +/- ', num2str(std(c_mc(1:100)))]) % Result for 100 iterations
fprintf('n')
disp('From MonteCarlo with 1000 iterations')
disp([num2str(mean(c_mc(1:1000))), ' +/- ', num2str(std(c_mc(1:1000)))]) % Result for 1000 iterations
fprintf('n')
disp('From MonteCarlo with 10000 iterations')
disp([num2str(mean(c_mc(1:10000))), ' +/- ', num2str(std(c_mc(1:10000)))]) % Result for 10000 iterations
fprintf('n')
disp('From MonteCarlo with 100000 iterations')
disp([num2str(mean(c_mc(1:100000))), ' +/- ', num2str(std(c_mc(1:100000)))]) % Result for 100000 iterations
fprintf('n')In this example, the uncertainty is propagated both analytically and with a Monte Carlo procedure, to then compare their results. First, the uncertainty arising on c is propagated analytically as the square root of the sum of the squared uncertainties of a and b.
Next, there's the Monte Carlo simulation. We randomly sample a and b from Gaussian distributions centered at their respective values and with standard deviations equal to their uncertainties. The samples values go into temporal varaibles a_mc and b_mc, whose subtraction is stored in the array c_mc. We repeat this process many times, here 100,000 times, and then compute the average and standard deviation on this array to estimate the central value adopted by c and its associated uncertainty. This example program reports the average and standard deviations after 10, 100, 1000, 10000 and 100000 iterations, from which we can evaluate how the estimates converge.
Here are the results of an example run:
Expected 4 +/- 2.2361
From Monte Carlo with 10 iterations 3.2695 +/- 3.1021
From Monte Carlo with 100 iterations 3.8985 +/- 2.5303
From Monte Carlo with 1000 iterations 3.9573 +/- 2.2683
From Monte Carlo with 10000 iterations 4.0182 +/- 2.2139
From Monte Carlo with 100000 iterations 3.9984 +/- 2.2439
As you see, both the result of the calculation and the error expected from analytical propagation (2.2361) are reasonably converged through iterations of the Monte Carlo procedure. Look especially at how the error, which is our main interest here, goes down and converges very close to the analytically propagated value. And check how even just 1000 iterations already produce a very reasonable estimate of the error.
Note that if you run the code yourself you will of course get different numbers, but the trend should always be there: with more iterations, your error should converge to values similar to those reported in the example (and to the analytically derived value).
Two more examples, now to estimate the uncertainty of parameters obtained by fitting data
In the previous section, we explored how Monte Carlo simulations can be used to propagate errors in simple arithmetic operations. Now, let's delve into more practical examples where we apply Monte Carlo simulations to estimate the uncertainties of parameters obtained by fitting data.
For this example, we will consider a linear fitting procedure. Given a set of data points (x, y), we want to find the best-fitting line of the form y = ax + b. However, real-world data often contains measurement errors or other uncertainties, which can impact the accuracy of the fitting parameters (a and b). To account for these uncertainties and obtain more reliable estimates, we can employ Monte Carlo simulations.
For this example we generate some sample data and add Gaussian noise centered around 0 with a width of +/- 1 to simulate measurement errors. This stands as our experimental vector of y values for the rest of the procedure.
clear; close all; clc
a = 5; b = 4; % Slope and intercept
x = [1:10]; % x from 1 to 10
y = a * x + b + normrnd(0, 1, 1, 10); % Simulate data and add Gaussian noise centered at 0 and with a width of +/- 1First, we perform linear regression using Matlab's fitlm function, which calculates the best-fit line and provides us with the analytically-derived slope and intercept along with their uncertainties:
analyticfit = fitlm(x,y) % fitlm performs a linear fit and displays the results, which include
% the analytically-derived slope and intercept and their uncertaintiesThe output of this procedure looks like this:
analyticfit =
Linear regression model:
Estimated Coefficients:
Estimate SE tStat pValue
________ _______ ______ __________
(Intercept) 3.6469 0.94449 3.8612 0.0048006
x1 5.0534 0.15222 33.198 7.3956e-10
Number of observations: 10, Error degrees of freedom: 8
Root Mean Squared Error: 1.38
R-squared: 0.993, Adjusted R-Squared: 0.992
F-statistic vs. constant model: 1.1e+03, p-value = 7.4e-10So, you see in this run and with the random numbers that the generator provided, the analytical fit yields slope = 5.0534 +/- 0.15222 and intercept = 3.6469 +/- 0.94449
Now, let's see what a Monte Carlo procedure finds. Just following the same idea of how we propagated the uncertainty of a subtraction operation, we must take the x and y vectors, add noise randomly (here only on y because we assume x has no error, but we could perfectly add noise to it too!), and we run a fitting procedure 10,000 times. This time we use polyfit, a Matlab function to fit polynomials (here of degree 1) that does not return uncertainties of the fitted parameters -which we'll get through Monte Carlo.
So, just like in the previous example, we add noise, run the linear fit, and store the parameters in vectors dedicated to the slopes (a_mc) and to the intercepts (b_mc):
a_mc = []; % Vector with slopes from Monte Carlo run
b_mc = []; % Vector with intercepts from Monte Carlo run
for i=1:10000 % 10k iterations, and in each:
yy = normrnd(y,1); % Create temporal vector where noise is added to the experimental y values
P = polyfit(x,yy,1); % Fit that vector to x using linear regresion as a polynomial of degree 1
b_mc = [b_mc; P(1)]; % Store all a values
a_mc = [a_mc; P(2)]; %
end
fprintf('n')
disp('From Monte Carlo with 10 iterations')
disp(['Intercept = ', num2str(mean(a_mc(1:10))), ' +/- ', num2str(std(a_mc(1:10)))])
disp(['Slope = ', num2str(mean(b_mc(1:10))), ' +/- ', num2str(std(b_mc(1:10)))])
fprintf('n')
disp('From Monte Carlo with 100 iterations')
disp(['Intercept = ', num2str(mean(a_mc(1:100))), ' +/- ', num2str(std(a_mc(1:100)))])
disp(['Slope = ', num2str(mean(b_mc(1:100))), ' +/- ', num2str(std(b_mc(1:100)))])
fprintf('n')
disp('From Monte Carlo with 1000 iterations')
disp(['Intercept = ', num2str(mean(a_mc(1:1000))), ' +/- ', num2str(std(a_mc(1:1000)))])
disp(['Slope = ', num2str(mean(b_mc(1:1000))), ' +/- ', num2str(std(b_mc(1:1000)))])
fprintf('n')
disp('From Monte Carlo with 10000 iterations')
disp(['Intercept = ', num2str(mean(a_mc(1:10000))), ' +/- ', num2str(std(a_mc(1:10000)))])
disp(['Slope = ', num2str(mean(b_mc(1:10000))), ' +/- ', num2str(std(b_mc(1:10000)))])
fprintf('n')Here's a set of example results I got:
From Monte Carlo with 10 iterations Intercept = 3.9435 +/- 0.97266 Slope = 5.0174 +/- 0.16053
From Monte Carlo with 100 iterations Intercept = 3.581 +/- 0.69941 Slope = 5.0716 +/- 0.11592
From Monte Carlo with 1000 iterations Intercept = 3.6461 +/- 0.69365 Slope = 5.0551 +/- 0.11206
From Monte Carlo with 10000 iterations Intercept = 3.6482 +/- 0.68509 Slope = 5.0533 +/- 0.10918
We observe that as the number of iterations increases, the Monte Carlo estimates for the intercept and slope converge towards the values obtained analytically with fitlm. The uncertainties obtained from Monte Carlo simulations (standard deviations) are consistent with the standard errors provided by fitlm, confirming the accuracy of the Monte Carlo method for propagating errors in the fitting parameters.
Yet, one of the main advantages of using a Monte Carlo method is that we can choose what the error in the inputs are, and we can not only consider uncertainties on the observations y but also on the dependent variables. These and other features of the method allow for a better, more controlled and knowledge-guided quantification and propagation of uncertainties through a calculation or algorithm. (See a dedicated section at the end that revolves around these and other advantages of the procedure, along with its disadvantages.)
An even harder example, from an actual application
The same approach we just used for linear procedures, can actually be applied to practically any kind of calculation, being particularly enticing for complex types of fitting procedures, where it can provide valuable insights into the reliability of model parameters and predictions.
While the above examples were mainly demonstrative, let me now show you one example coming from my own research where Monte Carlo simulations offer the only practical way to propagate errors into fitting parameters. The story begins with my need to fit a quite complex equation that looks like this:
The underlying data consists in an independent variable (our x) that is here (1/tcp), assumed of no noise, and a dependent variable (y) called R₂. My experimental dataset consists of multiple x, y values reported for different related systems that must all be fit together to produce a single number for kₑₓ, pA, pB, R₂A⁰, R₂B⁰ and d
