optimize.curve_fit不覆盖参数空间

问题描述

我正在尝试使用curve_fit将温度和降水数据拟合为周期函数。由于某些原因，curve_fit似乎无法测试bounds参数定义的整个参数空间。 我聚集了一个小模型来演示这一点。

############################################################
# First generate some data
import numpy as np

# Seed the random number generator for reproducibility
np.random.seed(0)
def test_func(x,b,c,d,e=0.0):
    return e*x + b * np.sin(2*np.pi * x/c + d)

def func_err(x,y):
    sumxy=0.0
    pos=0;
    for xtst in x:
        sumxy+=np.square(xtst-y[pos])
        pos+=1
    sumxy/=len(x)
    return sumxy

x_data = np.linspace(0,10,num=50)*2*np.pi
y_data = test_func(x_data,5.0,20,-0.20,0.1) + 1.0*np.random.normal(size=50)

# And plot it
import matplotlib.pyplot as plt

############################################################
# Now fit a simple sine function to the data
from scipy import optimize


params,params_covariance = optimize.curve_fit(test_func,x_data,y_data,p0=[1.0,18,0.0,0.0],bounds=([0.1,5,-5.0,-5.0],[100,100,5.0]))

print([params,func_err(y_data,test_func(x_data,params[0],params[1],params[2],params[3]))])

############################################################
# And plot the resulting curve on the data

plt.figure(figsize=(6,4))
plt.scatter(x_data,label='Data')
plt.plot(x_data,params[3]),label='Fitted function')

plt.legend(loc='best')

plt.show()

在给定p0=[1.0,0.0]的情况下，程序找到一个很好的解决方案，

但是使用p0=[1.0,0.0]之类的初始值，它会严重失败。

为什么例程无法覆盖边界给出的范围以找到其解决方案？

解决方法

我认为这是由于周期函数的性质所致。您的参数c确定函数的周期性。当您最初对周期性的猜测远离正确的周期性时，拟合将停留在局部最小值。

您可以认为，好像p0=[1.0,10,0.0,0.0]一样，拟合算法找到了局部最佳拟合，如第二个图所示，[ 0.65476428,11.14188385,-1.09652992,0.08971854]对[b,c,d,e]来说，它试图移动参数稍稍增加一点，但是周围的渐变表明这是最合适的，就像它位于参数空间的“谷”中一样，因此它在那里停止了迭代。

curve_fit不会探索您的整个参数空间：它仅从您最初的猜测开始（在这种情况下为p0，并使用启发式方法来查找局部最优值。

如果要浏览参数c的整个参数空间，可以实现简单的网格搜索。例如，您可以尝试在c的边界之间的所有值，并对每个curve_fit的值进行c，然后选择误差最小的一个。

这是示例代码：

def MSE(params,x_data,y_data):
    "to calclate mean square error,this is the same as your func_error"
    return ((test_func(x_data,*params)-y_data)**2).mean()

besterror = 10000
bestParam = None

for c_ in np.arange(5,100,1):
    # grid search for parameter c between 5 and 100,step size is 1.
    params,params_covariance = optimize.curve_fit(test_func,y_data,p0=[1.0,c_+0.5,0.0],bounds=([0.1,c_,-5.0,-5.0],[100,c_+1,5.0,5.0]))
    error = MSE(params,y_data)
    if error<besterror:
        besterror = error 
        bestParam = params

params = bestParam

print([params,func_err(y_data,test_func(x_data,params[0],params[1],params[2],params[3]))])

############################################################
# And plot the resulting curve on the data

params_covariance

plt.figure(figsize=(6,4))
plt.scatter(x_data,label='Data')
plt.plot(x_data,params[3]),label='Fitted function')

plt.legend(loc='best')

plt.show()

在这种情况下，无需在网格中搜索其他参数，因为curve_fit足以找到其他参数的最佳值。

这是一种蛮力方法，可能有一些库可以帮助您以更有效的方式进行操作。