Scipy.optimize.curve_fit 不符合余弦幂律
Scipy.optimize.curve_fit won't fit cosine power law
几个小时以来，我一直在尝试将模型拟合到（生成的）数据集，作为我一直在努力解决的问题的原因。我为函数 f(x) = A*cos^n(x)+b 生成了数据点，并添加了一些噪声。当我尝试使用此函数和 curve_fit 拟合数据集时，出现错误
./tester.py:10: RuntimeWarning: invalid value encountered in power
return Amp*(np.cos(x))**n + b
/usr/lib/python2.7/dist-packages/scipy/optimize/minpack.py:690: OptimizeWarning: Covariance of the parameters could not be estimated category=OptimizeWarning)
我用来生成数据点和拟合模型的代码如下：
#!/usr/bin/env python
from __future__ import print_function
import numpy as np
from scipy.optimize import curve_fit
from matplotlib.pyplot import figure, show, rc, plot
def f(x, Amp, n, b):
return np.real(Amp*(np.cos(x))**n + b)
x = np.arange(0, 6.28, 0.01)
randomPart = np.random.rand(len(x))-0.5
fig = figure()
sample = f(x, 5, 2, 5)+randomPart
frame = fig.add_subplot(1,1,1)
frame.plot(x, sample, label="Sample measurements")
popt, pcov = curve_fit(f, x, sample, p0=(1,1,1))
modeldata = f(x, popt[0], popt[1], popt[2])
print(modeldata)
frame.plot(x, modeldata, label="Best fit")
frame.legend()
frame.set_xlabel("x")
frame.set_ylabel("y")
show()
显示了噪声数据 - 请参见下图。
你们中有人知道发生了什么事吗？我怀疑它与进入复数域的幂律有关，因为函数的实部是 nowhere divergent。我试过只返回函数的实部，在 curve_fit 中设置实际边界，并使用 numpy 数组而不是 p0 的 python 列表。我是 运行 可用的 scipy 最新版本，scipy 0.17.0-1.
问题如下：
>>> (-2)**1.1
(-2.0386342710747223-0.6623924280875919j)
>>> np.array(-2)**1.1
__main__:1: RuntimeWarning: invalid value encountered in power
nan
与原生 python 浮点数不同，numpy 双精度数通常拒绝参与导致复杂结果的操作：
>>> np.sqrt(-1)
__main__:1: RuntimeWarning: invalid value encountered in sqrt
nan
作为快速解决方法，我建议向您的函数添加一个 np.abs 调用，并使用适当的边界进行拟合以确保这不会产生虚假拟合。如果您的模型接近真实并且您的样本（我的意思是样本中的余弦）为正，那么在它周围添加一个绝对值应该是一个空操作（更新：我意识到情况从来不是这样，请参阅正确的方法下面）。
def f(x, Amp, n, b):
return Amp*(np.abs(np.cos(x)))**n + b # only change here
有了这个小改动，我得到了这个：
供参考，拟合参数为 (4.96482314, 2.03690954, 5.03709923]) 与 (5,2,5) 相比。
经过深思熟虑后，我意识到余弦 总是 对于一半域为负 (duh)。所以我建议的解决方法可能有点问题，或者至少它的正确性很重要。另一方面，考虑包含 cos(x)^n 的原始公式，cos(x) 的负值只有在 n 是整数时才作为模型有意义，否则您会得到复杂的结果。既然无法解决Diophantine 拟合问题，就需要妥善处理。
最正确的方法（我的意思是最不可能使数据产生偏差的方法）是这样的：首先使用将数据转换为复数的模型进行拟合，然后在输出中采用复数：
def f(x, Amp, n, b):
return Amp*np.abs(np.cos(x.astype(np.complex128))**n) + b
这显然比我的解决方法效率低得多，因为在每个拟合步骤中我们都会创建一个新网格，并以复杂的算术和额外的幅度计算的形式做一些额外的工作。即使没有设置边界，这也给了我以下拟合：
参数为(5.02849409, 1.97655728, 4.96529108)。这些也很接近。但是，如果我们将这些值放回到实际模型中（没有 np.abs），我们会得到与 -0.37 一样大的虚部，这不是压倒性的但很重要。
所以第二步应该是用一个合适的模型重新拟合——一个有整数指数的模型。取你的拟合中明显的指数 2，并用这个模型做一个新的拟合。我不相信任何其他方法都能为您提供数学上合理的结果。也可以从原来的popt说起，希望确实接近真相。当然，我们可以使用一些柯里化的原始函数，但是使用模型的专用双特定版本要快得多。
from __future__ import print_function
import numpy as np
from scipy.optimize import curve_fit
from matplotlib.pyplot import subplots, show
def f_aux(x, Amp, n, b):
return Amp*np.abs(np.cos(x.astype(np.complex128))**n) + b
def f_real(x, Amp, n, b):
return Amp*np.cos(x)**n + b
x = np.arange(0, 2*np.pi, 0.01) # pi
randomPart = np.random.rand(len(x)) - 0.5
sample = f(x, 5, 2, 5) + randomPart
fig,(frame_aux,frame) = subplots(ncols=2)
for fr in frame_aux,frame:
fr.plot(x, sample, label="Sample measurements")
fr.legend()
fr.set_xlabel("x")
fr.set_ylabel("y")
# auxiliary fit for n value
popt_aux, pcov_aux = curve_fit(f_aux, x, sample, p0=(1,1,1))
modeldata = f(x, *popt_aux)
#print(modeldata)
print('Auxiliary fit parameters: {}'.format(popt_aux))
frame_aux.plot(x, modeldata, label="Auxiliary fit")
# check visually, test if it's close to an integer, but otherwise
n = np.round(popt_aux[1])
# actual fit with integral exponent
popt, pcov = curve_fit(lambda x,Amp,b,n=n: f_real(x,Amp,n,b), x, sample, p0=(popt_aux[0],popt_aux[2]))
modeldata = f(x, popt[0], n, popt[1])
#print(modeldata)
print('Final fit parameters: {}'.format([popt[0],n,popt[1]]))
frame.plot(x, modeldata, label="Best fit")
frame_aux.legend()
frame.legend()
show()
请注意，我更改了您的代码中的一些内容，但这并不真正影响我的观点。上面的图，所以既显示辅助配合又显示正确的：
输出：
Auxiliary fit parameters: [ 5.02628994 2.00886409 5.00652371]
Final fit parameters: [5.0288141074549699, 2.0, 5.0009730316739462]
重申一下：虽然辅助配合和正确配合之间可能没有视觉差异，但只有后者对您的问题给出了有意义的答案。
几个小时以来，我一直在尝试将模型拟合到（生成的）数据集，作为我一直在努力解决的问题的原因。我为函数 f(x) = A*cos^n(x)+b 生成了数据点，并添加了一些噪声。当我尝试使用此函数和 curve_fit 拟合数据集时，出现错误
./tester.py:10: RuntimeWarning: invalid value encountered in power
return Amp*(np.cos(x))**n + b
/usr/lib/python2.7/dist-packages/scipy/optimize/minpack.py:690: OptimizeWarning: Covariance of the parameters could not be estimated category=OptimizeWarning)
我用来生成数据点和拟合模型的代码如下：
#!/usr/bin/env python
from __future__ import print_function
import numpy as np
from scipy.optimize import curve_fit
from matplotlib.pyplot import figure, show, rc, plot
def f(x, Amp, n, b):
return np.real(Amp*(np.cos(x))**n + b)
x = np.arange(0, 6.28, 0.01)
randomPart = np.random.rand(len(x))-0.5
fig = figure()
sample = f(x, 5, 2, 5)+randomPart
frame = fig.add_subplot(1,1,1)
frame.plot(x, sample, label="Sample measurements")
popt, pcov = curve_fit(f, x, sample, p0=(1,1,1))
modeldata = f(x, popt[0], popt[1], popt[2])
print(modeldata)
frame.plot(x, modeldata, label="Best fit")
frame.legend()
frame.set_xlabel("x")
frame.set_ylabel("y")
show()
显示了噪声数据 - 请参见下图。
你们中有人知道发生了什么事吗？我怀疑它与进入复数域的幂律有关，因为函数的实部是 nowhere divergent。我试过只返回函数的实部，在 curve_fit 中设置实际边界，并使用 numpy 数组而不是 p0 的 python 列表。我是 运行 可用的 scipy 最新版本，scipy 0.17.0-1.
问题如下：
>>> (-2)**1.1
(-2.0386342710747223-0.6623924280875919j)
>>> np.array(-2)**1.1
__main__:1: RuntimeWarning: invalid value encountered in power
nan
与原生 python 浮点数不同，numpy 双精度数通常拒绝参与导致复杂结果的操作：
>>> np.sqrt(-1)
__main__:1: RuntimeWarning: invalid value encountered in sqrt
nan
作为快速解决方法，我建议向您的函数添加一个 np.abs 调用，并使用适当的边界进行拟合以确保这不会产生虚假拟合。如果您的模型接近真实并且您的样本（我的意思是样本中的余弦）为正，那么在它周围添加一个绝对值应该是一个空操作（更新：我意识到情况从来不是这样，请参阅正确的方法下面）。
def f(x, Amp, n, b):
return Amp*(np.abs(np.cos(x)))**n + b # only change here
有了这个小改动，我得到了这个：
供参考，拟合参数为 (4.96482314, 2.03690954, 5.03709923]) 与 (5,2,5) 相比。
经过深思熟虑后，我意识到余弦 总是 对于一半域为负 (duh)。所以我建议的解决方法可能有点问题，或者至少它的正确性很重要。另一方面，考虑包含 cos(x)^n 的原始公式，cos(x) 的负值只有在 n 是整数时才作为模型有意义，否则您会得到复杂的结果。既然无法解决Diophantine 拟合问题，就需要妥善处理。
最正确的方法（我的意思是最不可能使数据产生偏差的方法）是这样的：首先使用将数据转换为复数的模型进行拟合，然后在输出中采用复数：
def f(x, Amp, n, b):
return Amp*np.abs(np.cos(x.astype(np.complex128))**n) + b
这显然比我的解决方法效率低得多，因为在每个拟合步骤中我们都会创建一个新网格，并以复杂的算术和额外的幅度计算的形式做一些额外的工作。即使没有设置边界，这也给了我以下拟合：
参数为(5.02849409, 1.97655728, 4.96529108)。这些也很接近。但是，如果我们将这些值放回到实际模型中（没有 np.abs），我们会得到与 -0.37 一样大的虚部，这不是压倒性的但很重要。
所以第二步应该是用一个合适的模型重新拟合——一个有整数指数的模型。取你的拟合中明显的指数 2，并用这个模型做一个新的拟合。我不相信任何其他方法都能为您提供数学上合理的结果。也可以从原来的popt说起，希望确实接近真相。当然，我们可以使用一些柯里化的原始函数，但是使用模型的专用双特定版本要快得多。
from __future__ import print_function
import numpy as np
from scipy.optimize import curve_fit
from matplotlib.pyplot import subplots, show
def f_aux(x, Amp, n, b):
return Amp*np.abs(np.cos(x.astype(np.complex128))**n) + b
def f_real(x, Amp, n, b):
return Amp*np.cos(x)**n + b
x = np.arange(0, 2*np.pi, 0.01) # pi
randomPart = np.random.rand(len(x)) - 0.5
sample = f(x, 5, 2, 5) + randomPart
fig,(frame_aux,frame) = subplots(ncols=2)
for fr in frame_aux,frame:
fr.plot(x, sample, label="Sample measurements")
fr.legend()
fr.set_xlabel("x")
fr.set_ylabel("y")
# auxiliary fit for n value
popt_aux, pcov_aux = curve_fit(f_aux, x, sample, p0=(1,1,1))
modeldata = f(x, *popt_aux)
#print(modeldata)
print('Auxiliary fit parameters: {}'.format(popt_aux))
frame_aux.plot(x, modeldata, label="Auxiliary fit")
# check visually, test if it's close to an integer, but otherwise
n = np.round(popt_aux[1])
# actual fit with integral exponent
popt, pcov = curve_fit(lambda x,Amp,b,n=n: f_real(x,Amp,n,b), x, sample, p0=(popt_aux[0],popt_aux[2]))
modeldata = f(x, popt[0], n, popt[1])
#print(modeldata)
print('Final fit parameters: {}'.format([popt[0],n,popt[1]]))
frame.plot(x, modeldata, label="Best fit")
frame_aux.legend()
frame.legend()
show()
请注意，我更改了您的代码中的一些内容，但这并不真正影响我的观点。上面的图，所以既显示辅助配合又显示正确的：
输出：
Auxiliary fit parameters: [ 5.02628994 2.00886409 5.00652371]
Final fit parameters: [5.0288141074549699, 2.0, 5.0009730316739462]
重申一下：虽然辅助配合和正确配合之间可能没有视觉差异，但只有后者对您的问题给出了有意义的答案。