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

显示了嘈杂的数据 - 请参见下面的图片。

有人知道发生了什么吗？我怀疑这与幂律进入复数域有关，因为函数的实部无法发散。我已经尝试过仅返回函数的实部，设置曲线拟合的实际范围并使用numpy数组代替python列表p0。我正在运行最新版本的scipy，即scipy 0.17.0-1。

- Leander

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双精度浮点数通常不参与导致复杂结果的运算：

>>> 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)进行比较。

---

经过一番思考后，我意识到余弦函数对于你定义域的一半会永远是负数（呵呵）。因此，我建议的解决方法可能有些问题，或者至少它的正确性并不显然。另一方面，如果原公式包含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]

再次强调：虽然辅助适配和正确适配可能在视觉上没有区别，但只有后者才能给出有意义的答案解决您的问题。

- Andras Deak -- Слава Україні