Python中破幂律的曲线拟合答案

【问题标题】：Curve fitting with broken power law in PythonPython中破幂律的曲线拟合
【发布时间】：2015-11-23 02:49:27
【问题描述】：

我正在尝试遵循并重复使用由名为@ThePredator 的人建议的一段代码（使用我自己的数据）（我无法评论该线程，因为我目前没有所需的 50 声望）。完整代码如下：

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

我在运行这段代码时遇到了一些问题，我得到的错误如下：

y3 = (0)[x==xc]

TypeError: 'int' 对象没有属性 'getitem'

还有：

xc 未定义

我没有看到代码中缺少任何内容（不应该定义 xc 吗？）。

能否请作者 (@ThePredator) 或其他对此有所了解的人帮助我确定我没有看到的内容。

新版代码：

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

【问题讨论】：

如果你引用它的值，你确实需要定义xc（就像你在fitfunc中使用的比较一样。）y3 = (0)[x==xc]这一行不是有效的Python语法：@987654326 @ 是一个整数对象，不能使用 NumPy 样式的布尔索引进行索引，[x==xc]。
chi-sq 是无效的变量名，请改用下划线。
谢谢大家。理解了 y3 = (0)[x==xc]、chi-sq 和数组中的错误。关于变量xc，我应该用一个初始值来定义它，提供给拟合函数吗？
在答案的 cmets 中，您已经更改了代码，但没有在您的问题中更新它。您现在声称在生成的代码中有另一个未指定的错误。我还是不明白你是怎么从@ThePredator 明显的if(x==xc) y = 0; 变成y3 = (0)[x==xc] 的。
@ThePredator python 代码已经有y3 = (0)[x==xc] 它实际上是提问者提到if(x==xc) y = 0; 作为一段原始代码的原始问题。

标签： python-2.7 curve-fitting

【解决方案1】：

您的代码中有多个错误/错别字。

1) 您不能在 Python 中的变量名中使用 -（例如，卡方应该是 chi_square）

2) 您应该from numpy import array 或将array 替换为np.array。目前尚未定义名称array。

3) xc 没有定义，你应该在调用fitfunc()之前设置它。

4) y3 = (0)[x==xc] 无效，应该是（我认为）y3 = np.zeros(len(x))[x==xc] 或 y3 = np.zeros(np.sum(x==xc))

您对 fit_function() 的使用是错误的，因为它改变了图像的顺序。你想要的是：

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]

【讨论】：

感谢您的回复@M。马西亚斯。这是原始线程的链接：stackoverflow.com/questions/27153048/…
好的，现在它抱怨 y = np.concatenate((y1,y2,y3)) 的语法无效。
摆脱了所有的错误。 'y3 = np.zeros(np.sum(x==xc))' 行有效，而 'y3 = np.zeros(len(x)[x==xc]' 没有。所以代码现在运行没有错误，但它实际上并没有像预期的那样用破坏的幂律拟合数据。最佳拟合线是 y = 0 处的线。不知道如何修复它。有帮助吗？
好的，我会考虑一下，看看出了什么问题，或者我没有得到什么。如果您能在方便时提供帮助，那就太好了。
现在试试。似乎有问题。我正在尝试一些事情。将随着进度更新....

【解决方案2】：

嗨，请执行以下操作来定义您的功能，它会解决。 x 是一个数组（或列表），它应该将 y 作为一个数组（或列表）返回。然后你可以在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

【讨论】：