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

Wed Jun 8 10:33:45 EDT 2011

On Wed, Jun 8, 2011 at 10:27 AM, Johannes Radinger <JRadinger at gmx.at> wrote:
>
> -------- 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...

It doesn't really matter, if the scale only shows up in quadratic
terms, or as in my initial change I added a absolute value, so whether
it's positive or negative, it's still only one value, and we
interprete it as postive scale

s1 = sqrt(s1**2)

Josef

>
> 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
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>