## broken power law script

A forum for general discussion of the Python programming language.

### broken power law script

Hey evrybody... I am totally new to python... I found the following script... Please help me understand it...I know about MATLAB...either help me converting it to matlab or help me in understanding it so that i can run it...

Code: Select all
`import agpy.mpfit as mpfitimport numpy as npdef powerfit(xax,data,err=None,alphaguess=-2.0,scaleguess=1.0,quiet=True):    """    Fit a power law (a line in log-space) to data as a function of x    differs from 'plfit' because plfit fits a power law distribution,    this code simply fits a power law    """       logdata = np.log10(data)    if err is None: err = np.ones(data.shape,dtype='float')    def mpfitfun(data,err):        def f(p,fjac=None): return [0,np.ravel(((np.log10(p[0])+np.log10(xax)*p[1])-data)/err)]        return f           mp = mpfit.mpfit(mpfitfun(logdata,err),xall=[scaleguess,alphaguess],quiet=quiet)    fitp = mp.params    return fitp,mpdef brokenpowerfit(xax, data, err=None, alphaguess1=0.0, alphaguess2=-2.0, scaleguess=1.0,        breakpoint=None, quiet=True):    """    Fit a broken power law (a line in log-space) to data as a function of x    differs from 'plfit' because plfit fits a power law distribution,    this code simply fits a power law    This is a lot more intricate than the simple power law fit, since it involves    fitting two power laws with different slopes    Parameters:    p[0] - scale    p[1] - breakpoint    p[2] - power 1 (xax < breakpoint)    p[3] - power 2 (xax >= breakpoint)    There are 5 parameters (NOT 4) returned because there are two scales that are *NOT*    independent    returns: scale1,scale2,breakpoint,alpha1,alpha2    """       logdata = np.log10(data)    if err is None: err = np.ones(data.shape,dtype='float')    def brokenpowerlaw(p):        lowerhalf = (np.log10(p[0]) + np.log10(xax)*p[2]) * (xax < p[1])        # find the location at which both functions must be equal        scale2loc = np.argmin(np.abs(xax - p[1]))        scale2 = np.log10(xax[scale2loc])*(p[2] - p[3]) + np.log10(p[0])        upperhalf = (scale2 + np.log10(xax)*p[3]) * (xax >= p[1])        # DEBUG print "scale1: %15g   scale2: %15g  xaxind: %5i  xaxval: %15g  lower: %15g upper: %15g" % (p[0],scale2,scale2loc,np.log10(xax[scale2loc]),lowerhalf[scale2loc-1],upperhalf[scale2loc])        return lowerhalf+upperhalf    def mpfitfun(data,err):        def f(p,fjac=None): return [0,np.ravel((brokenpowerlaw(p)-data)/err)]        return f    if breakpoint is None: breakpoint = np.median(xax)    parinfo = [{}, {'mpminstep':xax.min(),'mpmaxstep':xax.max(),'step':xax.min()}, {}, {}]           mp = mpfit.mpfit(mpfitfun(logdata,err),xall=[scaleguess,breakpoint,alphaguess1,alphaguess2],quiet=quiet,parinfo=parinfo)    fitp = mp.params    scale2loc = np.argmin(np.abs(xax - fitp[1]))    scale2 = 10**( np.log10(xax[scale2loc])*(fitp[2] - fitp[3]) + np.log10(fitp[0]) )    fitp = np.array( [fitp[0],scale2] + fitp[1:].tolist() )    return fitp,mp`

I have a file with 3 columns(on which i have to run the script)... TIME MAGNITUDE MAGNITUDE_ERROR

how do we give the 3 columns as input in python...??

Thanks
Last edited by Yoriz on Sat Sep 14, 2013 11:29 pm, edited 1 time in total.

Posts: 1
Joined: Sat Sep 14, 2013 10:54 pm

### Re: broken power law script

What is the objective of your analysis? Do you only want to calculate the power fit? In MatLab you can do this via:

Code: Select all
`lobf = polyfit( log(x), log(y), 1 )a = exp( lobf(2) )b = lobf(1) % Assuming equation has form: y = a*x^b`

As far as inputting data in Python, you first must load the data in some usable way. e.g.

Code: Select all
`import numpy as npdata = open('textfile.txt','r').read().split('\n') # Have to remove header lines as appropriatedata = np.double( np.array( [r.split(delimiter) for r in data] ) )`

You now have an r x 3 array and can call the data as needed. e.g.

Code: Select all
`# Remember, Python operates on with zero as the first counter, not one like MatLabx1 = data[0,0]y1 = data[0,1]err = data[0,2]`
Python: 2.7 via Anaconda
Numpy: 1.7
Pandas: 0.11
OS: Windows 7
IDE: Spyder/IPython

tnknepp

Posts: 153
Joined: Mon Mar 11, 2013 7:41 pm