(* Uncertainty Propogation *) (* Author: Michael Ibrahim *) (* All of the code in this package is Copyright (C) December 1992 Michael Ibrahim 104 12th Ave. W. Virginia, MN 55792 wasfy@gac.edu and is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRENTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to Free Software Foundation, Inc. 675 Mass Ave Cambridge, MA 02139 USA. *) (* Version 1.01 *) (* Sources: `Data Reduction and Error Analysis for the Physical Sciences' by P. R. Bevington and D. K. Robinson. The algorithms for the symbolic and numerical propogation functions can be found in this book on pages 43 and 49-50 respectively. *) (* Warning: Redefines the $PreRead Command Comment out line 276 if this is not desired *) (* Examples: --------------------------------------- (* To propogate error simply wrap the function 'DoIt' around an expression containing uncertain data *) Uncertainty[1,2] - Uncertainty[2,3]//DoIt (* Uncertainty assumes that every datum is different unless the ID feature is used, for instance *) Uncertainty[1,2] + Uncertainty[1,2]//DoIt (* is evaluated assuming the two entries are different measured quantities. To over ride this feature simply use an ID, which can be any symbol, number, etc. *) Uncertainty[1,2,same] + Uncertainty[1,2,same]//DoIt (* In the following the first two quantities are actually identical, while the third is independent. *) Uncertainty[1,2,same] + Uncertainty[1,2,same] + Uncertainty[1,2]//DoIt (* If one wants to work with complicated expressions rather than retyping an ID simply let some variable be the quantity. Uncertainty.m only considers a quantity independent if it is the first time it is evaluated in Mathematica. *) a = Uncertainty[1,2]; b = Uncertainty[1,2]; a + a//DoIt a + b//DoIt (* The default propogation algorithm is a numerical method but a symbolic one is also included *) $PropogationFunction=SymbolicPropogation; a = Uncertainty[x,dx]; b = Uncertainty[y,dy]; c = a + b//DoIt (* The Uncertainty package also 'tracks' what each quantity depends on so it can warn you if you propogate 'dependent error' as 'independent error.' For instance, the following draws a Warning *) c + a//DoIt (* The correct quantity for the uncertainty can be found from *) a + a + b//DoIt (* A shorthand input form is also provided *) (2 +/- 3) / (4 +/- 5)//DoIt (* Uncertainty has the Listable attribute to simplify analysis with long lists of data. RelativeUncertainty assumes you enter the error as a percentage. It evaluates to head Uncertainty immediately. *) Uncertainty[Range[10],.5] - RelativeUncertainty[Range[10],10]//DoIt --------------------------------------- *) (* Making Different Propogation Functions: To make a propogation function simply complete the template below. "func" is a pure function whose arguments are all the quantities which have uncertainty. "args" is a list of the actual quantities of the function. "errors" is a list of the errors for the quantities. The arguments to "func" and the order of both "args" and "errors" are the same. One needs only to replace the '...' in the template below with code implementing the algorithm to find the uncertainty and placing it in the variable "err". PropFuncName[func_,args_,errors_,history_]:= Module[{val=func @@ args,err}, . . . Uncertainty[val,err,history] ] *) (* Bugs and Potential Improvements: It would be nice feature to have the Output form of 'Uncertainty' to only print out figures that are significant, e.g. (2.01 +/- .55) rather than (2.0124124 +/- .551123145). I should have the program check to see if the $PreRead variable is already defined before adding the +/- shorthand. I should also explore any weird affects due to precedence. NumericalPropogation function is lopsided. That is if you look up the algorithm you will see that it only checks the change in the function to the positive side of the data point. This could be bad if you have a function which is precipitous on the neg. side but flat on the other. A simple solution would be to check both sides and take the maximum error of the two. This should be trivial to do for whoever desires it. If anyone wants it done I could do it too. If one wraps 'DoIt' around an expression which has nothing with head Uncertainty a number of annoying error messages are made. *) BeginPackage["Uncertainty`Uncertainty`"] (**** USAGE MESSAGES ****) ValGrab::usage = "ValGrab[Uncertainty[]] obtains the \"value\" from an argument with head Uncertainty." ErrGrab::usage = "ErrGrab[Uncertainty[]] obtains the \"uncertainty\" from an argument with head Uncertainty." Uncertainty::usage = "Uncertainty[val,err] is a head given to an ordered pair of a value and the absolute uncertainty in that value. Uncertainty[val,err, ID] labels the quantity with ID so that it may be identified as unique from another quantity." RelativeUncertainty::usage = "RelativeUncertainty[val,err] is a head given to an ordered pair of a value and the percent uncertainty in that value. It immedediately evaluates to an expression with head Uncertainty." $PropogationFunction::usage = "$PropogationFunction is the algorithm which is used to propogate errors when the command DoIt is used. New propogation functions can be easily implemented. For details see the \"Making Different Propogation Functions\" blurb in the Uncertainty.m file." DoIt::usage = "DoIt[exp] performs the propogation of error on exp. Whenever an expression has the quantities with head \"Uncertainty\" simply follow the expression with \"//DoIt\" to evaluate." NumericalPropogation::usage = "NumericalPropogation is a function included in this package which performs numerical propogation of uncertainty. The algorithm can be found in Bevington. It is the default setting of $PropogationFunction." SymbolicPropogation::usage = "SymbolicPropogation is a function included in this package which performs numerical propogation of uncertainty. The algorithm can be found in Bevington. To use simply set $PropogationFunction=SymbolicPropogation." ID::usage = "ID is the head of the history of all quantities upon which an uncertain value was dependent. ID is used (as opposed to List) to prevent the Listable attribute of Uncertainty from threading over the history." IDGrab::usage = "IDGrab[Uncertainty[]] grabs the history (w/ head ID) from the quantity with head uncertainty." Begin["`Private`"] (**** AUXILIARY FUNCTIONS ****) func[a_,b_]:=Function[Evaluate[a],Evaluate[b]] (* functionify returns a list containing a pure function representing the expression, the arguments, the uncertainty in the arguments, and a history of dependents *) functionify[exp_]:= Module[{ vars, (* list of all Unique uncertainties *) dummies (* list of dummy variables for func *), comanc (* id's of anything that has common ancestry*), i (* just an index *)}, vars = Union[ (* Eliminates Duplicates *) Part[exp,Sequence @@ #]& /@ Position[exp,_Uncertainty]]; comanc = (Join @@ IDGrab /@ vars); (* all the history ID's in a list *) comanc = (Prepend[comanc, {}] //. (* all ID's appearing > 1 time *) ID[f_List,a___,x_,b___,x_,c___]:> (*Must be an easier way!*) ID[ Append[f,x],a,b,c ])[[1]]; Do[Message[ Uncertainty::comanc, Cases[vars, Uncertainty[_,_,ID[___,comanc[[i]],___]]] ], {i,Length[comanc]}]; dummies = Table[Unique[drk],{Length[vars]}]; Prepend[ MapAt[Union @@ #&, Transpose[{ValGrab[#],ErrGrab[#],IDGrab[#]}& /@ vars], 3], func[dummies,exp /. (Rule @@ #& /@ Transpose[{vars,dummies}])]]] (**** ERROR MESSAGES ****) Uncertainty::comanc = "Warning: It appears that the quantities `` have common ancestors. This may result in improper propogation of errors." Uncertainty::arg = "`` called with `` argument(s). 2 or 3 are expected." (**** EXPORTED FUNCTIONS ****) ValGrab[x_Uncertainty] := x[[1]] ErrGrab[x_Uncertainty] := x[[2]] IDGrab[x_Uncertainty] := x[[3]] Format[Uncertainty[a_,b_,_], OutputForm] := StringForm["(`` +/- ``)", a, b] Format[Uncertainty[a_,b_,_], TeXForm] := StringForm["{`` \pm ``}", a, b] Uncertainty[a_, b_] := Uncertainty[a,b,ID[Unique["Dubiety"]]] Uncertainty[a_, b_, c_] := Uncertainty[a,b,ID[c]] /; (Head[c] =!= ID) Uncertainty[a_] := Message[Uncertainty::arg, "Uncertainty", 1] Uncertainty[yo___] := Message[Uncertainty::arg, "Uncertainty", Length[{yo}]] /; Length[{yo}] =!= 3 RelativeUncertainty[a_,b_] := Uncertainty[a, a b/100] RelativeUncertainty[a_,b_,c_] := Uncertainty[a, a b/100, ID[c]] RelativeUncertainty[_] := Message[Uncertainty::arg, "RelativeUncertainty", 1] RelativeUncertainty[yo___] := Message[Uncertainty::arg, "RelativeUncertainty", Length[{yo}]] /; Length[{yo}] =!= 3 DoIt[exp_] := $PropogationFunction @@ functionify[exp] (** Different Propogation Algorithms -- Also Exported **) NumericalPropogation[func_,args_,errors_,history_]:= Module[{val=N[func @@ args]}, Uncertainty[val, N[Sqrt[Plus @@((Apply[func,#]& /@ (Plus[args,#]& /@DiagonalMatrix[errors])- val)^2) ]], history] ] SymbolicPropogation[fnc_,args_,errors_,history_]:= Module[{val = fnc @@ args, errs = errors, partials, err}, partials = #[fnc]& /@ (Derivative @@ #& /@ IdentityMatrix[Length[args]]); err = Sqrt[Plus @@ (((# @@ args & /@ partials) errs)^2)]; Uncertainty[val,err,history] ] $PropogationFunction=NumericalPropogation; (* Default Propogation Function *) (* This makes dealing with lists much more convenient *) SetAttributes[RelativeUncertainty,Listable] SetAttributes[Uncertainty,Listable] SetAttributes[DoIt,Listable] $PreRead = StringReplace[#, "+/-" -> "~Uncertainty~"]& (* Provides shorthand input form. Thanks J. Keiper and D. Withoff *) End[] Protect[Uncertainty] EndPackage[]