[SciPy-User] solving integration, density function

Warren Weckesser warren.weckesser at enthought.com
Thu Jan 6 07:10:14 EST 2011


On Thu, Jan 6, 2011 at 6:01 AM, Johannes Radinger <JRadinger at gmx.at> wrote:

> Thank you for the simplification of the formula,
> but I still get a different result in the case
> when x<s1 (eg. x=2, s1=3).
>
> here a code to try:
>
> ********************
> import math
> from scipy import stats
>
> s1 = 3
> m = 0
> p = 1
> x = 2
>
>

Your values are integers, so the division in the  expression (x-m)/s1 in the
code below uses integer division, and (x-m)/s1 will be 0.



> func = stats.norm.pdf(x, loc=m, scale=(s1))
> func2 = (1/(s1*math.sqrt(2*math.pi)) * math.exp(-0.5*((x-m)/s1)**2))
>
>

Change that to:

func2 = (1/(s1*math.sqrt(2*math.pi)) * math.exp(-0.5*((x-m)/float(s1))**2))

or change the behavior of the division operator by first executing

from __future__ import division


Warren



print func
> print func2
> ********************************
>
> /j
>
> -------- Original-Nachricht --------
> > Datum: Thu, 6 Jan 2011 05:48:25 -0600
> > Von: Warren Weckesser <warren.weckesser at enthought.com>
> > An: SciPy Users List <scipy-user at scipy.org>
> > Betreff: Re: [SciPy-User] solving integration, density function
>
> > On Thu, Jan 6, 2011 at 5:27 AM, Johannes Radinger <JRadinger at gmx.at>
> > wrote:
> >
> > > Hey
> > >
> > > Last time you helped me a lot with my normal
> > > probabilty density function. My problem now is
> > > quite simple, I think it's just a problem with
> > > the syntax (brackets):
> > >
> > > There are two ways to calculate the pdf, with the
> > > stats-function and with pure mathematically, but
> > > the give different results and I can't find the
> > > where I make the mistake:
> > >
> > >
> > > func1 = stats.norm.pdf(x, loc=m, scale=(s1))
> > > func2 =
> > > 1/((s1)*(math.sqrt(2*math.pi))))*(math.exp(((-0.5)*((x-m)/(s1)))**2)
> > >
> > > Where is the problem
> > >
> >
> >
> > func2 = 1/(s1*math.sqrt(2*math.pi)) * math.exp(-0.5*((x-m)/s1)**2)
> >
> >
> > Warren
> >
> >
> >
> > > thank you...
> > >
> > > /j
> > >
> > > -------- 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
> > > > _______________________________________________
> > > > SciPy-User mailing list
> > > > SciPy-User at scipy.org
> > > > http://mail.scipy.org/mailman/listinfo/scipy-user
> > >
> > > --
> > > NEU: FreePhone - kostenlos mobil telefonieren und surfen!
> > > Jetzt informieren: http://www.gmx.net/de/go/freephone
> > > _______________________________________________
> > > SciPy-User mailing list
> > > SciPy-User at scipy.org
> > > http://mail.scipy.org/mailman/listinfo/scipy-user
> > >
>
> --
> NEU: FreePhone - kostenlos mobil telefonieren und surfen!
> Jetzt informieren: http://www.gmx.net/de/go/freephone
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20110106/76b1e392/attachment.html>


More information about the SciPy-User mailing list