Do this rteally belong in sympy? What you do is using the classical delta-method approximation, today it is simple to get better results via simulation, like bootstrap.
Kjetil On Fri, Feb 8, 2013 at 3:01 PM, David Mashburn <[email protected]> wrote: > Hello sympy people! > > I couldn't find any information on the topic of "Error Analysis" in sympy > (possibly because terms like "standard", "error", and "propagation" don't > really catch this top well in search engines ;) > > Anyway, I just went ahead and wrote a basic converter from any expression to > the expression for the error; I just documented it in the docstring if you > are curious. > > It's a little rudimentary and could be better -- namely the error of (2*a) > is reported as 2*sqrt(ae**2/a**2)*Abs(a) instead of just 2*ae -- but it > seems to at least produce correct results. > > Please let me know if there is already a better way to do this or if you > find a bug with this code (or if you find the code useful)! > > Thanks! > -David Mashburn > > > def > GetCombinedErrorExpression(sympyExpression,variables,errorVariables,verbose=False,top=True): > '''Takes any symbolc expression with variables involving basic +-*/^ and > converts > it to the symbolic expression for the standard errors using error > propagation rules > > "variables" is the list of variables to use; anything else is assumed > to be constant > "errorVariables" is the corresponding list of variables for the > errors that will be used in the output > errorVariables are simple symbols used to generate the return > expression > verbose prints the expression and all sub-expressions generated by > recursion > > For cases not easily handlable, this simply returns None > > WARNING! No error handling is done in the case that there are > variables IN exponents; > you should only use error-free constants in exponents (but > bases can have errors) > > Simple examples: > > a,b,c,ae,be,ce = sympy.symbols('a,b,c,ae,be,ce') > > GetCombinedErrorExpression(a**2,[a],[ae]) > --> 2*a*ae > > GetCombinedErrorExpression(a+b+c,[a,b,c],[ae,be,ce]) ) > --> sqrt(ae**2 + be**2 + ce**2) > > GetCombinedErrorExpression(a*b*c,[a,b,c],[ae,be,ce]) ) > --> sqrt(ce**2/c**2 + be**2/b**2 + ae**2/a**2)*Abs(a*b*c) > > GetCombinedErrorExpression(a/b,[a,b],[ae,be]) > --> sqrt(be**2/b**2 + ae**2/a**2)*Abs(a/b) > > GetCombinedErrorExpression(a+2*b+c**2,[a,b,c],[ae,be,ce]) > (pretty form) --> > _____________________________ > / 2 2 > / 2 2 2 4*be *|b| > / ae + 4*c *ce + ---------- > / 2 > \/ b > ''' > if verbose: > print sympyExpression > if top: # Check that input variables are correct, but only do it once! > assert hasattr(variables,__len__) > assert hasattr(errorVariables.__len__) > assert len(variables)==len(errorVariables) > for v in variables: > assert v.__class__==sympy.Symbol > for v in errorVariables: > assert v.__class__==sympy.Symbol > > # Handle cases for sympyExpression (can be any of Add, Mul, Symbol, > Integer, Float, Number) > if sympyExpression.__class__==sympy.Add: > subExpressions = sympyExpression.as_coeff_add()[1] > subExpressionsError = [ > GetCombinedErrorExpression(s,variables,errorVariables,verbose,top=False) > for s in subExpressions ] > if None in subExpressionsError: > return None > else: > return sympy.sqrt(sum([s**2 for s in subExpressionsError])) > elif sympyExpression.__class__==sympy.Mul: > subExpressions = sympyExpression.as_coeff_mul()[1] > subExpressionsError = [ > GetCombinedErrorExpression(s,variables,errorVariables,verbose,top=False) > for s in subExpressions ] > if None in subExpressionsError: > return None > else: > errorsDivVars = [ subExpressionsError[i]/subExpressions[i] > for i in range(len(subExpressions)) ] > return sympy.Abs(sympyExpression)*sympy.sqrt(sum([s**2 for s in > errorsDivVars])) > elif sympyExpression.__class__==sympy.Pow: > subExpression,exponent = sympyExpression.as_base_exp() > subExpressionError = > GetCombinedErrorExpression(subExpression,variables,errorVariables,verbose,top=False) > if subExpressionError==None: > return None > return > sympyExpression*sympy.Abs(exponent)*subExpressionError/subExpression > elif sympyExpression.__class__ in (sympy.Symbol, > sympy.Integer, > sympy.Float, > sympy.Number): > if sympyExpression in variables: > return errorVariables[variables.index(sympyExpression)] # return > the appropriate error var > else: > return 0 # assume it's a constant > return sympyExpression > else: # kick any other cases as exceptions... > return None > > -- > You received this message because you are subscribed to the Google Groups > "sympy" group. > To unsubscribe from this group and stop receiving emails from it, send an > email to [email protected]. > To post to this group, send email to [email protected]. > Visit this group at http://groups.google.com/group/sympy?hl=en. > For more options, visit https://groups.google.com/groups/opt_out. > > -- "If you want a picture of the future - imagine a boot stamping on the human face - forever." George Orwell (1984) -- You received this message because you are subscribed to the Google Groups "sympy" group. To unsubscribe from this group and stop receiving emails from it, send an email to [email protected]. To post to this group, send email to [email protected]. Visit this group at http://groups.google.com/group/sympy?hl=en. For more options, visit https://groups.google.com/groups/opt_out.
