-1

Im trying to follow and re-use a piece of code (with my own data) suggested by someone named @ThePredator (I couldn't comment on that thread since I don't currently have the required reputation of 50). The full code is as follows:

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 = array([0.001, 0.524, 0.625, 0.670, 0.790, 0.910, 1.240, 1.640, 2.180, 35460])
y = array([7.435e-13, 3.374e-14, 1.953e-14, 3.848e-14, 4.510e-14, 5.702e-14, 5.176e-14, 6.0e-14,3.049e-14,1.12e-17])

""" 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()

Im having some problem running this code and the errors I get are as follows:

y3 = (0)[x==xc]

TypeError: 'int' object has no attribute 'getitem'

and also:

xc is undefined

I don't see anything missing in the code (xc shouldn't have to be defined?).

Could the author (@ThePredator) or someone else having knowledge about this please help me identify what i haven't seen.

  • New version of code:

    import numpy as np # This is the Numpy module
from scipy.optimize import curve_fit 
import matplotlib.pyplot as plt 

def fitfunc(x, a1, a2, xc):
    if x.all() < xc:
      y = x**a1
    elif x.all() > xc:
      y = x**(a1 - a2) * x**a2
    else:
      y = 0
    return y

xc = 2
x = np.array([0.001, 0.524, 0.625, 0.670, 0.790, 0.910, 1.240, 1.640, 2.180, 35460])
y = np.array([7.435e-13, 3.374e-14, 1.953e-14, 3.848e-14, 4.510e-14, 5.702e-14, 5.176e-14, 6.0e-14,3.049e-14,1.12e-17])

popt,pcov = curve_fit(fitfunc,x,y,p0=(1.0,1.0)) 

a1 = popt[0] 
a2 = popt[1]
residuals = y - fitfunc(x, a1, a2, xc)
chisq = sum((residuals**2)/fitfunc(x, a1, a2, xc)) 
curvey = [fitfunc(val, a1, a2, xc) for val in x] #  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()
Aqib
  • 11
  • 1
  • 4
  • You do need to define `xc` if you make reference to its value (as you do in the comparisons used in `fitfunc`.) The line `y3 = (0)[x==xc]` is not valid Python syntax: `(0)` is an integer object which cannot be indexed with NumPy-style boolean indexing, `[x==xc]`. – xnx Aug 28 '15 at 12:35
  • 1
    `chi-sq` is an invalid variable name, use underscore instead. – P. Camilleri Aug 28 '15 at 12:36
  • Thanks guys. Understood the mistakes in the line y3 = (0)[x==xc], chi-sq and the array. About the variable xc, should I define it with an initial value, which is supplied to the fitting function? – Aqib Aug 28 '15 at 12:53
  • 1
    In the comments to the answer you have changed your code, but haven't updated it in your question. You now claim to have another, unspecified error in the resulting code. I still can't figure out how you got from @ThePredator's obvious `if(x==xc) y = 0;` to `y3 = (0)[x==xc]`. – msw Aug 28 '15 at 13:15
  • @ThePredator python code already has `y3 = (0)[x==xc]` its actually the original question in which the asker mentioned `if(x==xc) y = 0;` as a piece of Origin code. – Aqib Aug 28 '15 at 13:19

2 Answers2

0

There are multiple errors/typos in your code.

1) You cannot use - in your variable names in Python (chi-square should be chi_square for example)

2) You should from numpy import array or replace array with np.array. Currently the name array is not defined.

3) xc is not defined, you should set it before calling fitfunc().

4) y3 = (0)[x==xc] is not valid, should be (I think) y3 = np.zeros(len(x))[x==xc] or y3 = np.zeros(np.sum(x==xc))

Your use of fit_function() is wrong, because it changes the order of the images. What you want is:

def fit_function(x, a1, a2, xc):
    if x < xc:
        y = x**a1
    elif x > xc:
        y = x**(a1 - a2) * x**a2
    else:
        y = 0
    return y
xc = 2 #or any value you want
curvey = [fit_function(val, a1, a2, xc) for val in x]
P. Camilleri
  • 12,664
  • 7
  • 41
  • 76
  • Thanks for your response @M. Massias. Heres the link to the original thread: http://stackoverflow.com/questions/27153048/implementing-a-broken-power-law-as-a-fitting-function-in-origin/ – Aqib Aug 28 '15 at 12:47
  • Ok now it complains about `y = np.concatenate((y1,y2,y3))` that its invalid syntax. – Aqib Aug 28 '15 at 12:58
  • Got rid of the all the mistakes. the line 'y3 = np.zeros(np.sum(x==xc))' worked whereas 'y3 = np.zeros(len(x)[x==xc]' did not. So the code is now running without errors, but its not actually fitting the data with a broken power law as expected. The best fit line is a line at y = 0. Not sure how to fix it. Any help? – Aqib Aug 28 '15 at 13:16
  • Ok I'll give it some thought to try and see whats going wrong or what is it that I'm not getting. It would be great if you could help when it is convenient for you. – Aqib Aug 28 '15 at 13:22
  • trying it now. There seems to be problems. I'm trying a few things. Will update with the progress.... – Aqib Aug 28 '15 at 14:08
  • @Aqib keep me posted about your progress – P. Camilleri Aug 31 '15 at 09:51
  • sorry I got drifted to something else. Got back to broken law now and Im still try to use the broken law function that you guided me to. Its not getting all four arguments that it required. Im editing my question to show the new code. Could you please comment what might be going wrong? – Aqib Sep 24 '15 at 14:42
  • you provided me with the function definition for the broken power law only, which I used with the code that I already have. I guess the code that I have, doesn't work properly. Im looking at it at the moment, but the issue is that Im a bit pressed on time. So if you could help, that'd be awesome. – Aqib Sep 28 '15 at 16:51
0

Hi Do the following to define your function, and it will solve. x is an array (or list) and it should return y as an array (or list). And then you can use it in curvefit.

def fit_function(x, a1, a2, xc):
    y = []
    for xx in x:
        if xx<xc:
            y.append(x**a1)
        elif xx>xc:
            y.append(x**(a1 - a2) * x**a2)
        else:
            y.append(0.0)
    return y
Vikas Chand
  • 31
  • 3