Confidence intervals in sympy.stats

[Buck Shlegeris]

I have the following code, using the statistics module. It answers questions like "If there's a 20% chance that my random variable is within 10 of the mean, what is the standard deviation?" with the call getStdDev(10,20).

def getStdDev(distance,confidence):
    x = Normal(0,1)
    confidenceDistance = x.confidence(confidence)[1]
    stddev = confidenceDistance * distance
    return stddev

The statistics module lets me do this with the confidence method. However, the statistics module is deprecated.

Is there an equivalent in the stats module?

[Matthew Rocklin]

This fails in solve

In [1]: from sympy.stats import *

In [2]: mu = Symbol('mu', real=True, bounded=True)

In [3]: sigma = Symbol('sigma', positive=True)

In [4]: X = Normal('X', mu, sigma)

In [5]: simplify(2 * P(X - mu < 10))

Out[5]: 

   ⎛    ___⎞    

   ⎜5⋅╲╱ 2 ⎟    

erf⎜───────⎟ + 1

   ⎝   σ   ⎠    

In [6]: solve(_ - 0.20, sigma)

NotImplementedError: 

No algorithms are implemented to solve equation -_Dummy_44 + erf(5*sqrt(2)/sigma)

Note that historically in statistics each distribution knew how to calculate these things on their own.  It was hard coded in.  In stats we express the problem with statistical operations and use SymPy to generate the right equations to solve in each case.  This is usually far more powerful and extensible but, in this particular case, performs badly.  The ideal path to solution here is to teach solve how to deal with error functions.

Does anyone know how hard this is?

--
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 sympy+un...@googlegroups.com.
To post to this group, send email to sy...@googlegroups.com.
Visit this group at http://groups.google.com/group/sympy.
For more options, visit https://groups.google.com/groups/opt_out.

[Aaron Meurer]

On Oct 6, 2013, at 10:13 PM, Matthew Rocklin <mroc...@gmail.com> wrote:

This fails in solve

In [1]: from sympy.stats import *

In [2]: mu = Symbol('mu', real=True, bounded=True)

In [3]: sigma = Symbol('sigma', positive=True)

In [4]: X = Normal('X', mu, sigma)

In [5]: simplify(2 * P(X - mu < 10))

Out[5]: 

   ⎛    ___⎞    

   ⎜5⋅╲╱ 2 ⎟    

erf⎜───────⎟ + 1

   ⎝   σ   ⎠    

In [6]: solve(_ - 0.20, sigma)

NotImplementedError: 

No algorithms are implemented to solve equation -_Dummy_44 + erf(5*sqrt(2)/sigma)

Note that historically in statistics each distribution knew how to calculate these things on their own.  It was hard coded in.  In stats we express the problem with statistical operations and use SymPy to generate the right equations to solve in each case.  This is usually far more powerful and extensible but, in this particular case, performs badly.  The ideal path to solution here is to teach solve how to deal with error functions.

Does anyone know how hard this is?

Does the error function have an inverse? The variable only appears once in the expression, so as long as each function is invertible, solving such expressions is easy. 

Aaron Meurer

[Matthew Rocklin]

Yes, erfinv

http://docs.sympy.org/dev/modules/functions/special.html#sympy.functions.special.error_functions.erfinv

[Aaron Meurer]

I guess solve() uses the inverse method, like

In [14]: tan(x).inverse()
Out[14]: sympy.functions.elementary.trigonometric.atan

So we just need to define that for erf and erfinv.

Aaron Meurer

[someone]

Hi,

> The ideal path to solution here is to teach solve
> how to deal with error functions.
>
> Does anyone know how hard this is?

I suggest to look at:

"A Simple Solver for Linear Equations Containing Nonlinear Operators"

by Norman Ramsey [1]. If done in right, this should work
for equations containing arbitrary (partially) invertible
special functions.

[1]: http://www.cs.tufts.edu/~nr/noweb/examples/solver.html