[SciPy-User] solving integration, density function

[Johannes Radinger]

-------- Original-Nachricht --------
> Datum: Tue, 21 Dec 2010 09:58:35 -0500
> Von: Skipper Seabold <jsseabold at gmail.com>
> An: SciPy Users List <scipy-user at scipy.org>
> Betreff: Re: [SciPy-User] solving integration, density function

> On Tue, Dec 21, 2010 at 9:33 AM, Johannes Radinger <JRadinger at gmx.at>
> wrote:
> >
> > -------- Original-Nachricht --------
> >> Datum: Tue, 21 Dec 2010 09:18:15 -0500
> >> Von: Skipper Seabold <jsseabold at gmail.com>
> >> An: SciPy Users List <scipy-user at scipy.org>
> >> Betreff: Re: [SciPy-User] solving integration, density function
> >
> >> On Tue, Dec 21, 2010 at 7:48 AM, Johannes Radinger <JRadinger at gmx.at>
> >> wrote:
> >> >
> >> > -------- Original-Nachricht --------
> >> >> Datum: Tue, 21 Dec 2010 13:20:47 +0100
> >> >> Von: Gregor Thalhammer <Gregor.Thalhammer at gmail.com>
> >> >> An: SciPy Users List <scipy-user at scipy.org>
> >> >> Betreff: Re: [SciPy-User] solving integration, density function
> >> >
> >> >>
> >> >> Am 21.12.2010 um 12:06 schrieb Johannes Radinger:
> >> >>
> >> >> > Hello,
> >> >> >
> >> >> > I am really new to python and Scipy.
> >> >> > I want to solve a integrated function with a python script
> >> >> > and I think Scipy should do that :)
> >> >> >
> >> >> > My task:
> >> >> >
> >> >> > I do have some variables (s, m, K,) which are now absolutely set,
> but
> >> in
> >> >> future I'll get the values via another process of pyhton.
> >> >> >
> >> >> > s = 400
> >> >> > m = 0
> >> >> > K = 1
> >> >> >
> >> >> > And have have following function:
> >> >> > (1/((s*K)*sqrt(2*pi)))*exp((-1/2*(((x-m)/s*K))^2) which is the
> >> density
> >> >> function of the normal distribution a symetrical curve with the mean
> >> (m) of
> >> >> 0.
> >> >> >
> >> >> > The total area under the curve is 1 (100%) which is for an
> >> integration
> >> >> from -inf to +inf.
> >> >> > I want to know x in the case of 99%: meaning that the integral (-x
> to
> >> >> +x) of the function is 0.99. Due to the symetry of the curve you can
> >> also set
> >> >> the integral from 0 to +x equal to (0.99/2):
> >> >> >
> >> >> > 0.99 =
> integral((1/((s*K)*sqrt(2*pi)))*exp((-1/2*(((x-m)/s*K))^2)),
> >> -x,
> >> >> x)
> >> >> > resp.
> >> >> > (0.99/2) =
> >> integral((1/((s*K)*sqrt(2*pi)))*exp((-1/2*(((x-m)/s*K))^2)),
> >> >> 0, x)
> >> >> >
> >> >> > How can I solve that question in Scipy/python
> >> >> > so that I get x in the end. I don't know how to write
> >> >> > the code...
> >> >>
> >> >>
> >> >> --->
> >> >> erf(x[, out])
> >> >>
> >> >>     y=erf(z) returns the error function of complex argument
> defined
> >> as
> >> >>     as 2/sqrt(pi)*integral(exp(-t**2),t=0..z)
> >> >> ---
> >> >>
> >> >> from scipy.special import erf, erfinv
> >> >> erfinv(0.99)*sqrt(2)
> >> >>
> >> >>
> >> >> Gregor
> >> >>
> >> >
> >> >
> >> > Thank you Gregor,
> >> > I only understand a part of your answer... I know that the integral
> of
> >> the density function is a error function and I know that the argument
> "from
> >> scipy.special import erf, erfinv" is to load the module.
> >> >
> >> > But how do I write the code including my orignial function so that I
> can
> >> modify it (I have also another function I want to integrate). how do i
> >> start? I want to save the whole code to a python-script I can then load
> e.g.
> >> into ArcGIS where I want to use the value of x for further
> calculations.
> >> >
> >>
> >> Are you always integrating densities?  If so, you don't want to use
> >> integrals probably, but you could use scipy.stats
> >>
> >> erfinv(.99)*np.sqrt(2)
> >> 2.5758293035489004
> >>
> >> from scipy import stats
> >>
> >> stats.norm.ppf(.995)
> >> 2.5758293035489004
> >
> >> Skipper
> >
> > The second function I want to integrate is different, it is a
> combination of two normal distributions like:
> >
> > 0.99 =
> integrate(0.6*((1/((s1*K)*sqrt(2*pi)))*exp((-1/2*(((x-m)/s1*K))^2))+0,4*((1/((s2*K)*sqrt(2*pi)))*exp((-1/2*(((x-m)/s2*K))^2)))
> >
> > and here again I know s1, s2, m and K and want to get x in the case when
> the integral is 0.99. What do I write into the script I want create?
> >
> > I think it is better if I can explain it with a graph but I don't know
> if I can just attach pictures to the mail-list-mail.
> >
> 
> The cdf is the integral of pdf and the ppf is the inverse of this
> function.  All of these functions can take an argument for loc and
> scale, which in your case loc=m and scale = s1*K.  I think you can get
> by with these.  You might be able to do something like this.  Not sure
> if this is correct with respect to your weightings, etc. I'd have to
> think more, but it might get you on the right track.
> 
> from scipy import optimize
> 
> def func(x,sigma1,sigma2,m):
>     return .6 * stats.norm.cdf(x, loc=m, scale=sigma1) + .4 *
> stats.norm.cdf(x,
>         loc=m, scale=sigma2) - .995
> 
> optimize.zeros.newton(func, 1., args=(s1*K,s2*K,m))
> 
ooh that is what I was looking for...you helped me a lot, I just need now to get the optimization subpackage to run...my IDE (IDLE always crashes when I want to run: "from scipy import optimize" ... I am working with Pyhton 2.6.5 and the actual Scipy package (downloaded today) on Windows 7...

but thanks a lot for all your help
/johannes



> Skipper
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user

-- 
GMX.at - Österreichs FreeMail-Dienst mit über 2 Mio Mitgliedern
E-Mail, SMS & mehr! Kostenlos: http://portal.gmx.net/de/go/atfreemail