Well, my opinion is that this sort of thing does belong in SymPy. Symbolic
solutions to these sorts of things can provide insight that simulations can
not or can provide intermediate representations to launch simulations in
the future.
SymPy can handle some of these manipulations through the stats module
In [1]: from sympy.stats import *
In [2]: mu = Symbol('mu', real=True)
In [3]: mu2 = Symbol('mu2', real=True)
In [4]: sigma = Symbol('sigma', positive=True)
In [5]: sigma2 = Symbol('sigma2', positive=True)
In [6]: x = Normal('X', mu, sigma)
In [7]: y = Normal('Y', mu2, sigma2)
In [8]: std(2*x)
Out[8]: 2⋅σ
In [9]: simplify(std(x + y))
Out[9]:
__________
╱ 2 2
╲╱ σ + σ₂
In [10]: simplify(std(y*x**2))
Out[10]:
__________________________________________________________
╱ 4 2 2 2 2 2 2 2 2 4 4 2
╲╱ μ ⋅σ₂ + 4⋅μ ⋅μ₂ ⋅σ + 6⋅μ ⋅σ ⋅σ₂ + 2⋅μ₂ ⋅σ + 3⋅σ ⋅σ₂
It supports other queries as well
In [18]: simplify(P(2*x > 1))
Out[18]:
⎛ ___ ⎞
⎜╲╱ 2 ⋅(2⋅μ - 1)⎟
erf⎜───────────────⎟
⎝ 4⋅σ ⎠ 1
──────────────────── + ─
2 2
On Fri, Feb 8, 2013 at 3:28 PM, Kjetil brinchmann Halvorsen <
[email protected]> wrote:
> 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.
>
>
>
--
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.
