使用 scipy.optmize.curve_fit 拟合日志数据

Question

我有两个变量x and y，我正在尝试使用来自scipy.optimize 的curve_fit 来拟合它们。

拟合数据的方程是 y=a(x^b) 形式的简单幂律。 I set the x and y axis to log scale（即ax.set_xscale('log') 和ax.set_yscale('log')）时的数据似乎很适合。

代码如下：

def fitfunc(x,p1,p2):
    y = p1*(x**p2)
    return y

popt_1,pcov_1 = curve_fit(fitfunc,x,y,p0=(1.0,1.0))
p1_1 = popt_1[0]
p1_2 = popt_1[1]
residuals1 = (ngal_mstar_1) - fitfunc(x,p1_1,p1_2)
xi_sq_1 = sum(residuals1**2) #The chi-square value

curve_y_1 = fitfunc(x,p1_1,p1_2) #This is the fit line seen in the graph

fig = plt.figure(figsize=(14,12))
ax1 = fig.add_subplot(111)
ax1.scatter(x,y,c='r')
ax1.plot(y,curve_y_1,'y.',linewidth=1)
ax1.legend(loc='best',shadow=True,scatterpoints=1)
ax1.set_xscale('log') #Scale is set to log
ax1.set_yscale('log') #SCale is set to log
plt.show()

当我对 x 和 y 使用真正的 log-log 值时，幂律拟合变为 y=10^(a+b*log(x))，即将右侧的幂提高到 10，因为它是 logbase 10。现在两者x 和 y 的值是 log(x) 和 log(y)。

上面的拟合似乎不太好。这是我使用的代码。

def fitfunc(x,p1,p2):
    y = 10**(p1+(p2*x))
    return y

popt_1,pcov_1 = curve_fit(fitfunc,np.log10(x),np.log10(y),p0=(1.0,1.0))

p1_1 = popt_1[0]
p1_2 = popt_1[1]
residuals1 = (y) - fitfunc((x),p1_1,p1_2)
xi_sq_1 = sum(residuals1**2)

curve_y_1 = fitfunc(np.log10(x),p1_1,p1_2) #The fit line uses log(x) here itself

fig = plt.figure(figsize=(14,12))
ax1 = fig.add_subplot(111)
ax1.scatter(np.log10(x),np.log10(y),c='r')
ax1.plot(np.log10(y),curve_y_1,'y.',linewidth=1)
plt.show()

两个图之间的唯一区别是拟合方程，而第二个图的值是独立记录的。我在这里做错了吗，因为我想要一个 log(x) vs log(y) 图和相应的拟合参数（斜率和截距）

Answer 1

您将幂律模型转换为对数对数是错误的，即您的第二次拟合实际上适合不同的模型。取你原来的模型y=a*(x^b)，两边取对数，你会得到log(y) = log(a) + b*log(x)。因此，您的对数尺度模型应该简单地读取y' = a' + b*x'，其中素数表示对数尺度的变量。该模型现在是一个线性函数，一个众所周知的结果是所有幂律都变成了 log-log 中的线性函数。

也就是说，您仍然可以期待您的两个版本之间存在一些小的差异，因为curve_fit 将优化最小二乘问题。因此，在对数尺度下，拟合将最小化拟合与数据之间的相对误差，而在线性尺度下，拟合将最小化绝对误差。因此，为了确定哪种方式实际上更适合您，您必须估计数据中的误差。您显示的数据当然在对数尺度上没有恒定的不确定性，因此在线性尺度上您的拟合可能更忠实。如果知道每个数据点中的错误的详细信息，那么您可以考虑使用sigma 参数。如果使用得当，这两种方法应该没有太大区别。在这种情况下，我更喜欢对数尺度拟合，因为模型更简单，因此可能在数值上更稳定。