[SciPy-User] solving integration, density function

Skipper Seabold jsseabold@gmail....
Tue Dec 21 08:58:35 CST 2010


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 :)
>> >> >
>> >> > 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))

Skipper


More information about the SciPy-User mailing list