Home › Code Snippets

Error propagation

Posted June 28th, 2010 by andyfaff
in
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 c = sin(A)
	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 c = cos(A)
	multithread dc = abs(-sin(A) * da)
End

Assuming you have

On February 14th, 2011 RGerkin says:

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

On February 14th, 2011 jjweimer says:

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

On February 14th, 2011 johnweeks says:

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

On February 14th, 2011 andyfaff says:

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.

Back to top