|
[Date Index]
[Thread Index]
[Author Index]
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
Prev by Date:
Re: Re: Mathematica 1
Next by Date:
Re: Mathematica 1
Previous by thread:
Mathematica on FreeBSD
Next by thread:
integer solution
|
