(* 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[]