## Error propagation

Average rating
(1 vote)

Some code for error propagation (use at own risk, not triple checked it).

```#pragma rtGlobals=1		// Use modern global access method.
//http://en.wikipedia.org/wiki/Propagation_of_uncertainty

Function EP_add(a, da, b, db, c, dc, [covar])
wave a, da, b, db, c, dc
variable covar
// C = A + B
duplicate/free a, tempa
duplicate/free b, tempb
duplicate/free da, tempda
duplicate/free db, tempdb

redimension/n=(dimsize(a, 0), dimsize(a, 1), dimsize(a, 2), dimsize(a, 3)) c, dc
if(paramisdefault(covar))
covar = 0
endif

multithread c = tempa + tempb
multithread dc = sqrt(tempda ^ 2 + tempdb ^ 2 + 2 * covar)
End

Function EP_sub(a, da, b, db, c, dc, [covar])
wave a, da, b, db, c, dc
variable covar
//C = A - B
duplicate/free a, tempa
duplicate/free b, tempb
duplicate/free da, tempda
duplicate/free db, tempdb

redimension/n=(dimsize(a, 0), dimsize(a, 1), dimsize(a, 2), dimsize(a, 3)) c, dc
if(paramisdefault(covar))
covar = 0
endif

multithread c = tempa - tempb
multithread dc = sqrt(tempda ^ 2 + tempdb ^ 2 - 2 * covar)
End

Function EP_mult(a, da, b, db, c, dc, [covar])
wave a, da, b, db, c, dc
variable covar
//C = A * B
duplicate/free a, tempa
duplicate/free b, tempb
duplicate/free da, tempda
duplicate/free db, tempdb

redimension/n=(dimsize(a, 0), dimsize(a, 1), dimsize(a, 2), dimsize(a, 3)) c, dc
if(paramisdefault(covar))
covar = 0
endif

multithread c = tempa * tempb
multithread dc = abs(c) * sqrt((tempda / tempa)^2 + (tempdb / tempb) ^ 2 + 2 * tempda * tempdb / (tempa * tempb) * covar)
End

Function EP_mulK(a, da, c, dc, k)
wave a, da, c, dc
variable k
//C = k * A
duplicate/free a, tempa
duplicate/free da, tempda

redimension/n=(dimsize(a, 0), dimsize(a, 1), dimsize(a, 2), dimsize(a, 3)) c, dc

multithread c = tempa * k
multithread dc = abs(tempda * k)
End

Function EP_div(a, da, b, db, c, dc, [covar])
wave a, da, b, db, c, dc
variable covar
//C =  A / B
duplicate/free a, tempa
duplicate/free b, tempb
duplicate/free da, tempda
duplicate/free db, tempdb

redimension/n=(dimsize(a, 0), dimsize(a, 1), dimsize(a, 2), dimsize(a, 3)) c, dc
if(paramisdefault(covar))
covar = 0
endif

multithread c = tempa / tempb
multithread dc = abs(c) * sqrt((tempda / tempa) ^ 2 + (tempdb / tempb) ^ 2 - 2 * tempda * tempdb / (tempa * tempb) * covar)
End

Function EP_pow(a, da, c, dc, k, [n])
wave a, da, c, dc
variable k, n
//C = n * (A ^ k)

duplicate/free a, tempa
duplicate/free da, tempda

redimension/n=(dimsize(a, 0), dimsize(a, 1), dimsize(a, 2), dimsize(a, 3)) c, dc
if(paramisdefault(n))
n = 1
endif

multithread c = n * (tempa ^ k)
multithread dc = abs(c * k * tempda / tempa)
End

Function EP_powK(a, da, c, dc, k, [n])
wave a, da, c, dc
variable k, n
//C = k ^ (A * n)
duplicate/free a, tempa
duplicate/free da, tempda

redimension/n=(dimsize(a, 0), dimsize(a, 1), dimsize(a, 2), dimsize(a, 3)) c, dc
if(paramisdefault(n))
n = 1
endif

multithread c = k ^ (tempA * n)
multithread dc = abs(c * n * tempda * ln(k))
End

Function EP_ln(a, da, c, dc, [k, n])
wave a, da, c, dc
variable k, n
//C = n * ln(k * A)
duplicate/free a, tempa
duplicate/free da, tempda

redimension/n=(dimsize(a, 0), dimsize(a, 1), dimsize(a, 2), dimsize(a, 3)) c, dc
if(paramisdefault(k))
k = 1
endif
if(paramisdefault(n))
n = 1
endif
multithread c = n * ln(k * tempa)
multithread dc = abs(n * tempda / tempa)
End

Function EP_log(a, da, c, dc, [k, n])
wave a, da, c, dc
variable k, n
//C = n * log(k * A)
duplicate/free a, tempa
duplicate/free da, tempda

redimension/n=(dimsize(a, 0), dimsize(a, 1), dimsize(a, 2), dimsize(a, 3)) c, dc
if(paramisdefault(k))
k = 1
endif
if(paramisdefault(n))
n = 1
endif
multithread c = n * log(k * tempa)
multithread dc = abs(n * tempda / (tempa * ln(10)))
End

Function EP_exp(a, da, c, dc, [k, n])
wave a, da, c, dc
variable k, n
//C = n * exp(k * A)
duplicate/free a, tempa
duplicate/free da, tempda

redimension/n=(dimsize(a, 0), dimsize(a, 1), dimsize(a, 2), dimsize(a, 3)) c, dc
if(paramisdefault(k))
k = 1
endif
if(paramisdefault(n))
n = 1
endif

multithread c = n * exp(k * tempa)
multithread dc = abs(k * tempda * c)
End

Function EP_sin(a, da, c, dc)
wave a, da, c, dc
//C = sin(A)
duplicate/free a, tempa
duplicate/free da, tempda

redimension/n=(dimsize(a, 0), dimsize(a, 1), dimsize(a, 2), dimsize(a, 3)) c, dc

multithread dc = abs(cos(A) * da)
End

Function EP_cos(a, da, c, dc)
wave a, da, c, dc
//C = cos(A)
duplicate/free a, tempa
duplicate/free da, tempda

redimension/n=(dimsize(a, 0), dimsize(a, 1), dimsize(a, 2), dimsize(a, 3)) c, dc

multithread dc = abs(-sin(A) * da)
End```

### Assuming you have

Assuming you have sufficiently fast and/or parallel computation, I wonder if it isn't usually easier to just get Monte Carlo error measurements. In other words, draw N initial values from a distribution that matches your measurements (e.g. mean = 3 cm, std = 0.4 cm), and run your usual, error-ignorant analysis on each of the N initial values, and then compute final errors from the distribution of the N results. Then you can use all of your existing code and never think about error propagation while you are writing new code.

### Isn't a principled drawback

Isn't a principled drawback to the MC method that it only gives the "true" values (those found by proper linear propagation methods) when N = infinity?

--
J. J. Weimer
Chemistry / Chemical & Materials Engineering, UAHuntsville

### Rick- check out Bootstrap

Rick- check out Bootstrap methods.

Statistical methods *all* have the drawback that they are exact only for N=infinity :)

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

### This code was developed for

This code was developed for reducing scattering data where the uncertainty in a detector pixel is treated as Poissonian, i.e. the square root of the count. Then one has to correct the data, dividing through by other detector images, integrating, subtracting, etc. There are a multitude of steps and at each point the uncertainty is propagated. However it takes about 4s to correct a single image. This would have to be repeated a bazillion of times (http://en.wikipedia.org/wiki/Indefinite_and_fictitious_numbers for MC errors), given that there are a huge number of steps during the reduction. This would be a long time . So I'm putting my faith in the Central Limit Theorem.