# [SciPy-User] solving integration, density function

josef.pktd@gmai... josef.pktd@gmai...
Tue Dec 21 09:50:54 CST 2010

```On Tue, Dec 21, 2010 at 9:58 AM, Skipper Seabold <jsseabold@gmail.com> wrote:
> On Tue, Dec 21, 2010 at 9:33 AM, Johannes Radinger <JRadinger@gmx.at> wrote:
>>
>> -------- Original-Nachricht --------
>>> Datum: Tue, 21 Dec 2010 09:18:15 -0500
>>> Von: Skipper Seabold <jsseabold@gmail.com>
>>> An: SciPy Users List <scipy-user@scipy.org>
>>> Betreff: Re: [SciPy-User] solving integration, density function
>>
>>> On Tue, Dec 21, 2010 at 7:48 AM, Johannes Radinger <JRadinger@gmx.at>
>>> wrote:
>>> >
>>> > -------- Original-Nachricht --------
>>> >> Datum: Tue, 21 Dec 2010 13:20:47 +0100
>>> >> Von: Gregor Thalhammer <Gregor.Thalhammer@gmail.com>
>>> >> An: SciPy Users List <scipy-user@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 :)
>>> >> >
>>> >> >
>>> >> > 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))

I think it's better to stick with the main namespace optimize.newton.
I think scipy.stats.distributions are using fsolve.

Skipper's way is the most direct way to calculate this.

Johannes,
scipy.stats.distribution has a lot of generic functions/methods to
work with distributions, and reading the source for some of them might
be useful to you.

Another alternative is to subclass stats.distribution and take
advantage of the generic methods directly. But compared to Skipper's
direct solution, this would be only worth it if you need more
properties. An old example of mine is at
http://mail.scipy.org/pipermail/scipy-user/2009-May/021182.html
(discussion was more for estimation)

Josef

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