sharing error propagation code with you (comments solicited)
- To: mathgroup at smc.vnet.net
- Subject: [mg61984] sharing error propagation code with you (comments solicited)
- From: Chris Chiasson <chris.chiasson at gmail.com>
- Date: Mon, 7 Nov 2005 02:10:41 -0500 (EST)
- Sender: owner-wri-mathgroup at wolfram.com
Dear MathGroup, I noticed that Interval and significance arithmetic tend to be pessimistic when it comes to propagation of error. I searched Math Source and found one package that deals with propagation of statistical error: http://library.wolfram.com/infocenter/MathSource/734/ That package has its own notation, and hard coded error propagation rules. I didn't test how it plays with mixed symbols and uncertain numbers though (which is important when working with units). It also can't handle Power with an uncertainty in the exponent combined with no uncertainty in the base, probably because the author(s) forgot that Power isn't Orderless like Times and Plus... Below, you will find code I wrote that automatically generates error propagation rules for functions, demonstrated with a few examples. It currently handles Exp, Log, Sin, Cos, Tan, Sec, Csc, Erf, Plus, Times, Subtract, and Power. Extending it to another function is a matter of adding the function to a list. It plays well with units. I have also included a comparison of its error propagation with Interval and significance arithmetic. What modifications should I make? { Cell[CellGroupData[{ Cell[BoxData[ RowBox[{"If", "[", RowBox[{ RowBox[{"errorstuffisset", "===", "True"}], ",", "Null", ",", RowBox[{ RowBox[{"Block", "[", RowBox[{ RowBox[{"{", RowBox[{"g", ",", RowBox[{"singleparamfuns", "=", RowBox[{"{", RowBox[{ "Exp", ",", "Log", ",", "Sin", ",", "Cos", ",", "Tan", ",", "Sec", ",", "Csc", ",", "Erf"}], "}"}]}], ",", RowBox[{"multiparamfuns", "=", RowBox[{"{", RowBox[{ "Plus", ",", "Times", ",", "Subtract", ",", "Power"}], "}"}]}]}], "}"}], ",", RowBox[{ RowBox[{"FullSimplify", "[", RowBox[{"{", RowBox[{ RowBox[{ RowBox[{"g", "[", RowBox[{"error", ",", RowBox[{"errorophead", "[", RowBox[{ RowBox[{"error", "[", RowBox[{"x1_", ",", "dx1_"}], "]"}], ",", RowBox[{"error", "[", RowBox[{"x2_", ",", "dx2_"}], "]"}]}], "]"}], ",", RowBox[{ RowBox[{"error", "[", RowBox[{ RowBox[{"errorophead", "[", RowBox[{"xr1", ",", "xr2"}], "]"}], ",", RowBox[{"Sqrt", "[", RowBox[{ SuperscriptBox[ RowBox[{"(", RowBox[{ RowBox[{"D", "[", RowBox[{ RowBox[{"errorophead", "[", RowBox[{"xr1", ",", "xr2"}], "]"}], ",", "xr1"}], "]"}], " ", "dxr1"}], ")"}], "2"], "+", SuperscriptBox[ RowBox[{"(", RowBox[{ RowBox[{"D", "[", RowBox[{ RowBox[{"errorophead", "[", RowBox[{"xr1", ",", "xr2"}], "]"}], ",", "xr2"}], "]"}], " ", "dxr2"}], ")"}], "2"]}], "]"}]}], "]"}], "/.", RowBox[{"{", RowBox[{ RowBox[{"xr1", "\[Rule]", "x1"}], ",", RowBox[{"dxr1", "\[Rule]", "dx1"}], ",", RowBox[{"xr2", "\[Rule]", "x2"}], ",", RowBox[{ "dxr2", "\[Rule]", "dx2"}]}], "}"}]}]}], "]"}], "/.", RowBox[{"(", RowBox[{ RowBox[{ RowBox[{"{", RowBox[{ "errorophead", "\[Rule]", "#"}], "}"}], "&"}], "/@", "multiparamfuns"}], ")"}]}], ",", RowBox[{ RowBox[{"g", "[", RowBox[{"error", ",", RowBox[{"errorophead", "[", RowBox[{ RowBox[{"error", "[", RowBox[{"x1_", ",", "dx1_"}], "]"}], ",", RowBox[{"x2_", "?", "NumericQ"}]}], "]"}], ",", RowBox[{ RowBox[{"error", "[", RowBox[{ RowBox[{"errorophead", "[", RowBox[{"xr1", ",", "xr2"}], "]"}], ",", RowBox[{"Sqrt", "[", RowBox[{ SuperscriptBox[ RowBox[{"(", RowBox[{ RowBox[{"D", "[", RowBox[{ RowBox[{"errorophead", "[", RowBox[{"xr1", ",", "xr2"}], "]"}], ",", "xr1"}], "]"}], " ", "dxr1"}], ")"}], "2"], "+", SuperscriptBox[ RowBox[{"(", RowBox[{ RowBox[{"D", "[", RowBox[{ RowBox[{"errorophead", "[", RowBox[{"xr1", ",", "xr2"}], "]"}], ",", "xr2"}], "]"}], " ", "dxr2"}], ")"}], "2"]}], "]"}]}], "]"}], "/.", RowBox[{"{", RowBox[{ RowBox[{"xr1", "\[Rule]", "x1"}], ",", RowBox[{"dxr1", "\[Rule]", "dx1"}], ",", RowBox[{"xr2", "\[Rule]", "x2"}], ",", RowBox[{"dxr2", "\[Rule]", "0"}]}], "}"}]}]}], "]"}], "/.", RowBox[{"(", RowBox[{ RowBox[{ RowBox[{"{", RowBox[{ "errorophead", "\[Rule]", "#"}], "}"}], "&"}], "/@", "multiparamfuns"}], ")"}]}], ",", RowBox[{ RowBox[{"g", "[", RowBox[{"error", ",", RowBox[{"errorophead", "[", RowBox[{ RowBox[{"x1_", "?", "NumericQ"}], ",", RowBox[{"error", "[", RowBox[{"x2_", ",", "dx2_"}], "]"}]}], "]"}], ",", RowBox[{ RowBox[{"error", "[", RowBox[{ RowBox[{"errorophead", "[", RowBox[{"xr1", ",", "xr2"}], "]"}], ",", RowBox[{"Sqrt", "[", RowBox[{ SuperscriptBox[ RowBox[{"(", RowBox[{ RowBox[{"D", "[", RowBox[{ RowBox[{"errorophead", "[", RowBox[{"xr1", ",", "xr2"}], "]"}], ",", "xr1"}], "]"}], " ", "dxr1"}], ")"}], "2"], "+", SuperscriptBox[ RowBox[{"(", RowBox[{ RowBox[{"D", "[", RowBox[{ RowBox[{"errorophead", "[", RowBox[{"xr1", ",", "xr2"}], "]"}], ",", "xr2"}], "]"}], " ", "dxr2"}], ")"}], "2"]}], "]"}]}], "]"}], "/.", RowBox[{"{", RowBox[{ RowBox[{"xr1", "\[Rule]", "x1"}], ",", RowBox[{"dxr1", "\[Rule]", "0"}], ",", RowBox[{"xr2", "\[Rule]", "x2"}], ",", RowBox[{ "dxr2", "\[Rule]", "dx2"}]}], "}"}]}]}], "]"}], "/.", RowBox[{"(", RowBox[{ RowBox[{ RowBox[{"{", RowBox[{ "errorophead", "\[Rule]", "#"}], "}"}], "&"}], "/@", RowBox[{"{", RowBox[{ "Plus", ",", "Times", ",", "Subtract", ",", "Power"}], "}"}]}], ")"}]}], ",", "\[IndentingNewLine]", RowBox[{ RowBox[{"g", "[", RowBox[{"error", ",", RowBox[{"errorophead", "[", RowBox[{"error", "[", RowBox[{"x1_", ",", "dx1_"}], "]"}], "]"}], ",", RowBox[{ RowBox[{"error", "[", RowBox[{ RowBox[{ "errorophead", "[", "xr1", "]"}], ",", RowBox[{"Sqrt", "[", SuperscriptBox[ RowBox[{"(", RowBox[{ RowBox[{"D", "[", RowBox[{ RowBox[{ "errorophead", "[", "xr1", "]"}], ",", "xr1"}], "]"}], " ", "dxr1"}], ")"}], "2"], "]"}]}], "]"}], "/.", RowBox[{"{", RowBox[{ RowBox[{"xr1", "\[Rule]", "x1"}], ",", RowBox[{ "dxr1", "\[Rule]", "dx1"}]}], "}"}]}]}], "]"}], "/.", RowBox[{"(", RowBox[{ RowBox[{ RowBox[{"{", RowBox[{ "errorophead", "\[Rule]", "#"}], "}"}], "&"}], "/@", "singleparamfuns"}], ")"}]}]}], "}"}], "]"}], "/.", RowBox[{"g", "\[Rule]", "TagSetDelayed"}]}]}], "]"}], ";", RowBox[{"errorstuffisset", "=", "True"}], ";"}]}], "]"}]], "Input", CellLabel->"In[1]:="], Cell[BoxData[ RowBox[{ RowBox[{"General", "::", "\<\"spell\"\>"}], RowBox[{ ":", " "}], "\<\"Possible spelling error: new symbol name \ \\\"\\!\\(dxr1\\)\\\" is similar to existing symbols \\!\\({dx1, xr1}\ \\). \\!\\(\\*ButtonBox[\\\"More\[Ellipsis]\\\", \ ButtonStyle->\\\"RefGuideLinkText\\\", ButtonFrame->None, \ ButtonData:>\\\"General::spell\\\"]\\)\"\>"}]], "Message", CellLabel->"From In[1]:="], Cell[BoxData[ RowBox[{ RowBox[{"General", "::", "\<\"spell\"\>"}], RowBox[{ ":", " "}], "\<\"Possible spelling error: new symbol name \ \\\"\\!\\(dxr2\\)\\\" is similar to existing symbols \\!\\({dx2, xr2}\ \\). \\!\\(\\*ButtonBox[\\\"More\[Ellipsis]\\\", \ ButtonStyle->\\\"RefGuideLinkText\\\", ButtonFrame->None, \ ButtonData:>\\\"General::spell\\\"]\\)\"\>"}]], "Message", CellLabel->"From In[1]:="] }, Open ]], Cell[BoxData[{ RowBox[{ RowBox[{ RowBox[{"reps", "[", "0", "]"}], "=", RowBox[{"error", "\[Rule]", "PlusMinus"}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{ RowBox[{"reps", "[", "1", "]"}], "=", RowBox[{"x_", ":>", RowBox[{"PlusMinus", "[", RowBox[{"x", ",", RowBox[{"x", " ", SuperscriptBox["10", RowBox[{"-", RowBox[{"Precision", "[", "x", "]"}]}]]}]}], "]"}]}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{ RowBox[{"reps", "[", "2", "]"}], "=", RowBox[{ RowBox[{"Interval", "[", "x_", "]"}], "\[RuleDelayed]", RowBox[{"PlusMinus", "[", RowBox[{ RowBox[{"Mean", "[", "x", "]"}], ",", RowBox[{ RowBox[{"Last", "[", "x", "]"}], "-", RowBox[{"Mean", "[", "x", "]"}]}]}], "]"}]}]}], ";"}]}], "Input", CellLabel->"In[102]:="], Cell[CellGroupData[{ Cell[BoxData[{ RowBox[{ RowBox[{ RowBox[{"error", "[", RowBox[{"100", ",", "0.01"}], "]"}], "+", RowBox[{"error", "[", RowBox[{"200", ",", "0.01"}], "]"}]}], "/.", RowBox[{"reps", "[", "0", "]"}]}], "\[IndentingNewLine]", RowBox[{ RowBox[{ RowBox[{"SetPrecision", "[", RowBox[{"100", ",", RowBox[{"-", RowBox[{"Log", "[", RowBox[{"10", ",", RowBox[{"0.01", "/", "100"}]}], "]"}]}]}], "]"}], "+", RowBox[{"SetPrecision", "[", RowBox[{"200", ",", RowBox[{"-", RowBox[{"Log", "[", RowBox[{"10", ",", RowBox[{"0.01", "/", "200"}]}], "]"}]}]}], "]"}]}], "/.", RowBox[{"reps", "[", "1", "]"}]}], "\[IndentingNewLine]", RowBox[{ RowBox[{ RowBox[{"Interval", "[", RowBox[{"100", "+", RowBox[{"{", RowBox[{ RowBox[{"-", "0.01"}], ",", "0.01"}], "}"}]}], "]"}], "+", RowBox[{"Interval", "[", RowBox[{"200", "+", RowBox[{"{", RowBox[{ RowBox[{"-", "0.01"}], ",", "0.01"}], "}"}]}], "]"}]}], "/.", RowBox[{"reps", "[", "2", "]"}]}]}], "Input", CellLabel->"In[105]:="], Cell[BoxData[ RowBox[{ "300", "\[PlusMinus]", "0.01414213562373095`"}]], "Output", CellLabel->"Out[105]="], Cell[BoxData[ RowBox[{ "300.`4.176091259055681", "\[PlusMinus]", "0.019999999999999997`"}]], "Output", CellLabel->"Out[106]="], Cell[BoxData[ RowBox[{ "300.`", "\[PlusMinus]", "0.020000000000095497`"}]], "Output", CellLabel->"Out[107]="] }, Open ]], Cell[CellGroupData[{ Cell[BoxData[ RowBox[{ RowBox[{"a", " ", RowBox[{"Erf", "[", RowBox[{"1", "/", SqrtBox[ RowBox[{"error", "[", RowBox[{"2", ",", "0.0001"}], "]"}]]}], "]"}]}], "/.", RowBox[{"reps", "[", "0", "]"}]}]], "Input", CellLabel->"In[125]:="], Cell[BoxData[ RowBox[{"a", " ", RowBox[{"(", RowBox[{ RowBox[{"Erf", "[", FractionBox["1", SqrtBox["2"]], "]"}], "\[PlusMinus]", "0.000012098536225957168`"}], ")"}]}]], "Output", CellLabel->"Out[125]="] }, Open ]], Cell[CellGroupData[{ Cell[BoxData[ RowBox[{ RowBox[{"(*", RowBox[{"Family", " ", "Count"}], "*)"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"Fold", "[", RowBox[{"ReplaceAll", ",", RowBox[{"You", "+", "Spouse", "+", "Kids", "+", "Pets"}], ",", RowBox[{"{", RowBox[{ RowBox[{"{", RowBox[{ RowBox[{"You", "\[Rule]", "1"}], ",", RowBox[{"Spouse", "\[Rule]", "1"}], ",", RowBox[{"Kids", "\[Rule]", RowBox[{"error", "[", RowBox[{ RowBox[{"2", "+", FractionBox["1", "2"]}], ",", FractionBox["1", "2"]}], "]"}]}], ",", RowBox[{"Pets", "\[Rule]", RowBox[{"error", "[", RowBox[{"1", ",", "1"}], "]"}]}]}], "}"}], ",", RowBox[{"reps", "[", "0", "]"}]}], "}"}]}], "]"}], "\[IndentingNewLine]", RowBox[{"%", "//", "N"}]}]}]], "Input", CellLabel->"In[122]:="], Cell[BoxData[ RowBox[{ FractionBox["11", "2"], "\[PlusMinus]", FractionBox[ SqrtBox["5"], "2"]}]], "Output", CellLabel->"Out[122]="], Cell[BoxData[ RowBox[{ "5.5`", "\[PlusMinus]", "1.118033988749895`"}]], "Output", CellLabel->"Out[123]="] }, Open ]] } -- http://chrischiasson.com/contact/chris_chiasson