Scipy.optimize.curve_fit不符合余弦幂律

Question

几个小时以来,我一直试图将模型拟合到(生成的)数据集中,作为我一直在努力解决的问题的一个问题.我为函数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()

显示噪声数据 - 请参见下图.

你们当中有没有人知道发生了什么？我怀疑它与进入复杂领域的幂律有关,因为函数的真实部分无处不同.我已经尝试仅返回函数的实部,在curve_fit中设置实际边界并使用numpy数组而不是p0的python列表.我正在运行最新版本的scipy,scipy 0.17.0-1.

Answer 1

问题如下:

>>> (-2)**1.1
(-2.0386342710747223-0.6623924280875919j)
>>> np.array(-2)**1.1
__main__:1: RuntimeWarning: invalid value encountered in power
nan

与原生的python浮标不同,numpy double通常拒绝参与操作,导致复杂的结果:

>>> 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是整数,否则您将得到复杂的结果.由于我们无法解决丢番图配件问题,我们需要妥善处理.

最合适的方式(我指的是最不可能偏向数据的方式)是这样的:首先使用将数据转换为复数的模型进行拟合,然后在输出上采用复杂的幅度:

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]

重申一下:虽然辅助配合与正确配合之间可能没有视觉差异,但只有后者才能为您的问题提供有意义的答案.