[SciPy-User] How to fit a curve/function?

Wed Jun 8 10:27:43 EDT 2011

-------- Original-Nachricht --------
> Datum: Wed, 8 Jun 2011 10:12:58 -0400
> Von: josef.pktd at gmail.com
> An: SciPy Users List <scipy-user at scipy.org>
> Betreff: Re: [SciPy-User] How to fit a curve/function?

> On Wed, Jun 8, 2011 at 9:41 AM, Johannes Radinger <JRadinger at gmx.at>
> wrote:
> >
> > -------- Original-Nachricht --------
> >> Datum: Wed, 8 Jun 2011 07:10:38 -0400
> >> Von: josef.pktd at gmail.com
> >> An: SciPy Users List <scipy-user at scipy.org>
> >> Betreff: Re: [SciPy-User] How to fit a curve/function?
> >
> >> On Wed, Jun 8, 2011 at 6:52 AM, Johannes Radinger <JRadinger at gmx.at>
> >> wrote:
> >> > Hello,
> >> >
> >> > I've got following function describing any kind of animal dispersal
> >> kernel:
> >> >
> >> > def pdf(x,s1,s2):
> >> >    return
> >>
> (p/(math.sqrt(2*math.pi*s1**2))*numpy.exp(-((x-0)**(2)/(2*s1**(2)))))+((1-p)/(s2*math.sqrt(2*math.pi))*numpy.exp(-((x-0)**(2)/(2*s2**(2)))))
> >> >
> >> > On the other hand I've got data from literature with which I want to
> fit
> >> the function so that I get s1, s2 and x.
> >> > Ususally the data in the literature are as follows:
> >> >
> >> > Example 1: 50% of the animals are between -270m and +270m and 90%
>  are
> >> between -500m and + 500m
> >> >
> >> > Example 2: 84% is between - 5000 m and +5000m, and 73% are between
> >> -1000m and +1000m
> >> >
> >> > So far as I understand an integration of the function is needed to
> solve
> >> for s1 and s2 as all the literature data give percentage (area under
> the
> >> curve) Can that be used to fit the curve or can that create ranges for
> s1
> >> and s2.
> >>
> >> I don't see a way around integration.
> >>
> >> If you have exactly 2 probabilities, then you can you a solver like
> >> scipy.optimize.fsolve to match the probabilites
> >> eg.
> >> 0.5 = integral pdf from -270 to 270
> >> 0.9 = integral pdf from -500 to 500
> >>
> >> If you have more than 2 probabilities, then using optimization of a
> >> weighted function of the moment conditions would be better.
> >>
> >> Josef
> >
> >
> >
> > Hello again
> >
> > I tried following, but without success so far. What do I have to do
> excactly...
> >
> > import numpy
> > from scipy import stats
> > from scipy import integrate
> > from scipy.optimize import fsolve
> > import math
> >
> > p=0.3
> >
> > def pdf(x,s1,s2):
> >    return
> (p/(math.sqrt(2*math.pi*s1**2))*numpy.exp(-((x-0)**(2)/(2*s1**(2)))))+((1-p)/(s2*math.sqrt(2*math.pi))*numpy.exp(-((x-0)**(2)/(2*s2**(2)))))
> >
> > def equ(s1,s2):
> >    0.5==integrate.quad(pdf,-270,270,args=(s1,s2))
> >    0.9==integrate.quad(pdf,-500,500,args=(s1,s2))
> >
> > result=fsolve(equ, 1,500)
> >
> > print result
> 
> equ needs to return the deviation of the equations (I changed some
> details for s1 just to try it)
> 
> import numpy
> from scipy import stats
> from scipy import integrate
> from scipy.optimize import fsolve
> import math
> 
> p=0.3
> 
> def pdf(x,s1,s2):
>     return
> (p/(math.sqrt(2*math.pi*s1**2))*numpy.exp(-((x-0)**(2)/(2*s1**(2)))))+((1-p)/(math.sqrt(2*math.pi*s2**2))*numpy.exp(-((x-0)**(2)/(2*s2**(2)))))
> 
> def equ(arg):
>     s1,s2 = numpy.abs(arg)
>     cond1 = 0.5 - integrate.quad(pdf,-270,270,args=(s1,s2))[0]
>     cond2 = 0.9 - integrate.quad(pdf,-500,500,args=(s1,s2))[0]
>     return [cond1, cond2]
> 
> result=fsolve(equ, [200., 1200])
> 
> print result
> 
> but in the results I get the parameters are very close to each other
> [-356.5283675   353.82544075]
> 
> the pdf looks just like a mixture of 2 normals both with loc=0, then
> maybe the cdf of norm can be used directly

Thank you for that hint... First yes these are 2 superimposed normals but for other reasons I want to use the original formula instead of the stats.functions...

anyway there is still a thing...the locator s1 and s2 are like the scale parameter of stats.norm so the are both + and -. For fsolve above it seems that I get only one parameter (s1 or s2) but for the positive and negative side of the distribution. So in actually there are four parameters -s1, +s1, -s2, +s2. How can I solve that? Maybe I can restrict the fsolve to look for the two values only in the positive range...

any guesses?

/J

> 
> >>> from scipy import stats
> >>> stats.norm.cdf(270, scale=350) - stats.norm.cdf(-270, scale=350)
> 0.55954705470577526
> >>>
> >>> stats.norm.cdf(270, scale=354) - stats.norm.cdf(-270, scale=354)
> 0.55436474670960978
> >>> stats.norm.cdf(500, scale=354) - stats.norm.cdf(-500, scale=354)
> 0.84217642881921018
> 
> Josef
> >
> >
> > /Johannes
> >>
> >> >
> >> > /Johannes
> >> >
> >> > --
> >> > NEU: FreePhone - kostenlos mobil telefonieren!
> >> > 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
> >> >
> >> _______________________________________________
> >> SciPy-User mailing list
> >> SciPy-User at scipy.org
> >> http://mail.scipy.org/mailman/listinfo/scipy-user
> >
> > --
> > NEU: FreePhone - kostenlos mobil telefonieren!
> > 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
> >
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user

-- 
NEU: FreePhone - kostenlos mobil telefonieren!			
Jetzt informieren: http://www.gmx.net/de/go/freephone