Monte Carlo Error Propagation

Error propagation method for an arbitrary analytic function with different error types

Error propagation is of central interest in modern science and in most cases done by assuming gaussian errors for the parameters and the calculating the partial derivatives (see https://en.wikipedia.org/wiki/Propagation_of_uncertainty#Simplification).
This method, however, has three major drawbacks:
a) it is only exactly true for linear functions or functions that can well be approximated by a linear function, but breaks down completely for example in case of f(a,b) = a/b when the ratio becomes small, while the error remains significant (see example 2).
b) in the simple version it is impossible to combine parameters which have different error distributions that a gaussian distribution (e.g. binomial)
c) for complex functions the calculation of partial derivatives can be tedious
We wrote a simple Monte Carlo based error propagation, which allows to prevent all of these drawbacks. The file example.m contains a variety of different concrete examples how to use the method and shows where the simple Gaussian error propagation method breaks down. Thanks to modern computers this method allows an exact error propagation by numerical Monte Carlo parameter generation.
The method essentially consists of two functions: generateMCparameters and propagateErrorWithMC
The first part generates an distribution of MC parameter values with the following options:
% errorType: gaussian, binomial, bootstrapMean, bootstrapDistribution
% params: depends on errorType (gauss: x,dx; binomial: n,k; bootstrap: array measured values
% plot: (optional) plot final distribution
% numSamples: (optional) number of MC samples

at the current state the following four distributions can be generated:
- Gaussian: specified by the mean value and the sigma
- Binomial: defined by n and k
- bootstrapMean: this was implemented, because a lot of times one measures a signal and knows it has a mean value, but the readings fluctuate, in this case the user can enter the measured values (x_1,x_2,....x_n) and using bootstrapping a distribution centered around the mean is generated (see also https://en.wikipedia.org/wiki/Bootstrapping_(statistics)#Estimating_the_distribution_of_sample_mean).
- bootstrapDistribution: in case the measured values themselves fluctuate (and not the just the reading) this method allows to generate a MC array directly from the measured values (x_1,x_2,....x_n).

once the distributions of the parameters are generated one can propagate them. For example in case of the function f(a,b) = a/b:
A = generateMCparameters('gaussian',[2,0.2]);
B = generateMCparameters('gaussian',[0.5,0.2]);
paramMatrix = [A;B];
funToProp = @(x) x(1)./x(2);
[funValue,funCI,funSamples] = propagateErrorWithMC(funToProp, paramMatrix);

with the following options:
% funOfInterest function that should be evaluated
% params: matrix of column vectors, each row represents sampled parameters
% CIthreshold: (optional) confidence interval threshold, default: 0.68
% plot: (optional) plot final distribution
% method: (optional) method to determine funValue (median (default), mean, maximum)

The default value for the confidence interval is CIthreshold = 0.68. The CI is then determined by integrating the function value distribution from +/- inf until the value reaches (1-CIthreshold)/2.

The final plot shows the values within the CI in green and everything outside in blue, also a exponential fit is performed to compare the final distribution with a gaussian.

==== Version 1.0 (2016-07-14) ====
upload of the initial version written by:
Carsten Robens and Stefan Brakhane

known bugs/needs improvement:
- The CI from propagateErrorWithMC for a purely binomial distribution with small n does not reproduce the well known Clopper Pearson CI (see https://en.wikipedia.org/wiki/Binomial_proportion_confidence_interval#Clopper-Pearson_interval).