Implementing a Broken Power Law as a fitting funct

Question

Good day, I'm trying to use the function builder in origin (OriginLab) to create a new function to fit with, the broken power law (http://en.wikipedia.org/wiki/Power_law#Broken_power_law)

So, I think I got the actual function part down. For this I used

if(x<xc)
    y =x^a1;
if(x>xc)
    y = x^(a1-a2)*x^a2;
if(x==xc)
    y = 0;

Where xc, a1 and a2 are the parameters. However, I then get to the point where you have to choose a bunch of stuff (parameter ranges, script you run to guess initial values, etc) and I have no clue what to put in there. Does anyone have some experience with this?

Answer 1

Even though the question asks for a suggestion using OriginLab, this one is using Python as the OP has accepted to give it a try!

The curve-fitting method that exists in Python is from the Scipy package (curve_fit). All python packages for windows can be downloaded in a jiffy from THIS WEBSITE HERE!

When doing a curve fit, the first thing one needs to know is the fitting equation. As you already know the broken power law equation that will fit your data, you are ready to go.

The code for fitting an example data, lets call them x and y. Here the fitting parameters are a1 and a2.

import numpy as np # This is the Numpy module
from scipy.optimize import curve_fit # The module that contains the curve_fit routine
import matplotlib.pyplot as plt # This is the matplotlib module which we use for plotting the result

""" Below is the function that returns the final y according to the conditions """
def fitfunc(x,a1,a2):
    y1 = (x**(a1) )[x<xc]
    y2 = (x**(a1-a2) )[x>xc]
    y3 = (0)[x==xc]
    y = np.concatenate((y1,y2,y3))
    return y

x = Your x data here
y = Your y data here

""" In the above code, we have imported 3 modules, namely Numpy, Scipy and matplotlib """

popt,pcov = curve_fit(fitfunc,x,y,p0=(10.0,1.0)) #here we provide random initial parameters a1,a2

a1 = popt[0] 
a2 = popt[1]
residuals = y - fitfunc(x,a1,a2)
chi-sq = sum( (residuals**2)/fitfunc(x,a1,a2) ) # This is the chi-square for your fitted curve

""" Now if you need to plot, perform the code below """
curvey = fitfunc(x,a1,a2) # This is your y axis fit-line

plt.plot(x, curvey, 'red', label='The best-fit line')
plt.scatter(x,y, c='b',label='The data points')
plt.legend(loc='best')
plt.show()

Just insert your data here and it should work just fine!! For more detailed information if needed on how the code works, CHECK OUT THIS WEBSITE I couldn't find a proper example for your fit-function, so x and y have been left blank. But once you have your data, just plug them in!