如何设置curve_fit的初始值以寻找最佳优化而不仅仅是局部优化？

Question

如何设置curve_fit的初始值以寻找最佳优化而不仅仅是局部优化？

pythonscikit-learnscipyleast-squaresbest-fit-curve

3

我正在尝试拟合一个幂律函数，以找到最佳拟合参数。但是，如果初始参数猜测不同，则“最佳拟合”输出也会不同。除非我找到正确的初始猜测，否则我只能得到局部优化而不是最佳优化。有没有办法找到**适当的初始猜测**？下面是我的代码。请随意提出任何建议。谢谢！

import numpy as np
import pandas as pd
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
%matplotlib inline

# power law function
def func_powerlaw(x,a,b,c):
    return a*(x**b)+c

test_X = [1.0,2,3,4,5,6,7,8,9,10]
test_Y =[3.0,1.5,1.2222222222222223,1.125,1.08,1.0555555555555556,1.0408163265306123,1.03125, 1.0246913580246915,1.02]

predict_Y = []
for x in test_X:
    predict_Y.append(2*x**-2+1)

如果我使用默认的初始猜测值，即 p0 = [1,1,1] 进行对齐

popt, pcov = curve_fit(func_powerlaw, test_X[1:], test_Y[1:], maxfev=2000)


plt.figure(figsize=(10, 5))
plt.plot(test_X, func_powerlaw(test_X, *popt),'r',linewidth=4, label='fit: a=%.4f, b=%.4f, c=%.4f' % tuple(popt))
plt.plot(test_X[1:], test_Y[1:], '--bo')
plt.plot(test_X[1:], predict_Y[1:], '-b')
plt.legend()
plt.show()

拟合结果如下，不是最佳拟合。

如果我将初始猜测值更改为p0 = [0.5,0.5,0.5]

popt, pcov = curve_fit(func_powerlaw, test_X[1:], test_Y[1:], p0=np.asarray([0.5,0.5,0.5]), maxfev=2000)

我可以得到最佳拟合

---------------------2018年10月7日更新-------------------------------------------------------------------------------------------------------------------------

由于我需要运行数千甚至数百万个Power Law函数，使用@James Phillips的方法太昂贵了。除了curve_fit之外，什么方法是适当的？例如sklearn，np.linalg.lstsq等。

- Zed Fang

3个回答

3

并没有简单的答案：如果有的话，它就会被实现在 curve_fit 中，那么它就不必要求你提供起始点了。一个合理的方法是先拟合齐次模型 y = a*x**b。假设 y 为正（这通常是使用幂律时的情况），可以通过粗略而快速的方法来完成此操作：在对数-对数尺度上，log(y) = log(a) + b*log(x)，这是线性回归，可以用 np.linalg.lstsq 解决。这给出了 log(a) 和 b 的候选值；使用此方法的 c 的候选值为 0。

test_X = np.array([1.0,2,3,4,5,6,7,8,9,10])
test_Y = np.array([3.0,1.5,1.2222222222222223,1.125,1.08,1.0555555555555556,1.0408163265306123,1.03125, 1.0246913580246915,1.02])

rough_fit = np.linalg.lstsq(np.stack((np.ones_like(test_X), np.log(test_X)), axis=1), np.log(test_Y))[0]
p0 = [np.exp(rough_fit[0]), rough_fit[1], 0]

结果是您在第二张图片中看到的很好的匹配。

顺便说一句，最好一次将test_X制成NumPy数组。否则，您会首先切片X[1:]，这将作为整数数组转换为NumPy阵列，然后使用负指数抛出错误。(我想1.0的目的是将其制成浮点数数组？这就是应该使用dtype=np.float参数的原因。)

- user6655984

谢谢您的回复。在这种情况下，我需要适配超过数千甚至数百万个幂律函数。我认为首先应该应用对数-对数比例尺是一个好的解决方案，我也在其他地方找到了这个解决方案。我只是想知道，在使用对数-对数比例尺后，您是否认为使用 sklearn 包中的线性模型 可以解决这个问题？ - Zed Fang

3

除了Stack Overflow的Welcome和James Phillips所提供的非常好的答案之外，即“没有易于普遍适用的方法”和“差分进化通常有助于找到良好的起始点（甚至是良好的解决方案！）如果比curve_fit()慢一些”，让我单独回答一下，您可能会发现这很有帮助。

首先，curve_fit()默认任何参数值都是一个令人沮丧的坏主意。这种行为没有任何可能的理由，您和其他人应该将参数具有默认值视为curve_fit()实现中的严重错误，并假装不存在此错误。永远不要相信这些默认值是合理的。

从数据的简单绘图中，显然a=1，b=1，c=1是非常糟糕的起始值。函数衰减，因此b < 0。事实上，如果您从a=1，b=-1，c=1开始，就会找到正确的解决方案。

它可能有助于对参数设置合理的边界。即使将 c 的边界设置为（-100，100）也可能有所帮助。与 b 的符号一样，我认为您可以从数据的简单绘图中看到该边界。当我尝试在您的问题上进行此操作时，如果初始值为 b = 1 ，则对 c 的边界不起作用，但对于 b = 0 或 b = -5 。

更重要的是，虽然您在图中打印了最佳拟合参数 popt ，但您没有打印存储在 pcov 中的变量之间的不确定性或相关性，因此您对结果的解释是不完整的。如果您查看了这些值，您会发现从 b = 1 开始不仅会导致糟糕的值，而且还会导致参数的巨大不确定性和非常高的相关性。这是拟合告诉您找到了一个糟糕的解决方案。不幸的是，从curve_fit返回的pcov很难解包。

让我推荐一下 lmfit (https://lmfit.github.io/lmfit-py/)（免责声明：我是主要开发者之一）。除了其他功能外，该模块强制你提供非默认初始值，并更轻松地生成一个更完整的报告。对于您的问题，即使从 a=1, b=1, c=1 开始，也可以提供更有意义的指示表明存在问题：

from lmfit import Model
mod = Model(func_powerlaw)
params = mod.make_params(a=1, b=1, c=1)
ret = mod.fit(test_Y[1:], params, x=test_X[1:])
print(ret.fit_report())

这将打印出：

[[Model]]
    Model(func_powerlaw)
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 1318
    # data points      = 9
    # variables        = 3
    chi-square         = 0.03300395
    reduced chi-square = 0.00550066
    Akaike info crit   = -44.4751740
    Bayesian info crit = -43.8835003
[[Variables]]
    a: -1319.16780 +/- 6892109.87 (522458.92%) (init = 1)
    b:  2.0034e-04 +/- 1.04592341 (522076.12%) (init = 1)
    c:  1320.73359 +/- 6892110.20 (521839.55%) (init = 1)
[[Correlations]] (unreported correlations are < 0.100)
    C(a, c) = -1.000
    C(b, c) = -1.000
    C(a, b) =  1.000

这是 a = -1.3e3 +/- 6.8e6 -- 定义不太清楚！此外，所有参数都完全相关。

将 b 的初始值更改为-0.5：

params = mod.make_params(a=1, b=-0.5, c=1) ## Note !
ret = mod.fit(test_Y[1:], params, x=test_X[1:])
print(ret.fit_report())

给予

[[Model]]
    Model(func_powerlaw)
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 31
    # data points      = 9
    # variables        = 3
    chi-square         = 4.9304e-32
    reduced chi-square = 8.2173e-33
    Akaike info crit   = -662.560782
    Bayesian info crit = -661.969108
[[Variables]]
    a:  2.00000000 +/- 1.5579e-15 (0.00%) (init = 1)
    b: -2.00000000 +/- 1.1989e-15 (0.00%) (init = -0.5)
    c:  1.00000000 +/- 8.2926e-17 (0.00%) (init = 1)
[[Correlations]] (unreported correlations are < 0.100)
    C(a, b) = -0.964
    C(b, c) = -0.880
    C(a, c) =  0.769

这略微更好。

简而言之，初始值总是很重要的，结果不仅包括最佳匹配值，还包括不确定性和相关性。

- M Newville

1

谢谢您的回复！这是一个很好的提醒。我之前没有真正看过pcov。我有一个与此相关的快速问题。我尝试了初始猜测[0.5，0.5，0.5]和[1.0，-1.0，1.0]，两者都给出了优化点。然而，pcov的结果略有不同~~您知道为什么吗？ - Zed Fang

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- James Phillips · Accepted Answer

以下是使用scipy.optimize.differential_evolution遗传算法的示例代码，涉及您的数据和方程。该scipy模块使用拉丁超立方体算法来确保对参数空间进行彻底搜索，因此需要搜索范围，本示例中这些范围基于数据最大和最小值。对于其他问题，如果您知道参数值的范围，则可能需要提供不同的搜索范围。

import numpy, scipy, matplotlib
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy.optimize import differential_evolution
import warnings

# power law function
def func_power_law(x,a,b,c):
    return a*(x**b)+c

test_X = [1.0,2,3,4,5,6,7,8,9,10]
test_Y =[3.0,1.5,1.2222222222222223,1.125,1.08,1.0555555555555556,1.0408163265306123,1.03125, 1.0246913580246915,1.02]


# function for genetic algorithm to minimize (sum of squared error)
def sumOfSquaredError(parameterTuple):
    warnings.filterwarnings("ignore") # do not print warnings by genetic algorithm
    val = func_power_law(test_X, *parameterTuple)
    return numpy.sum((test_Y - val) ** 2.0)


def generate_Initial_Parameters():
    # min and max used for bounds
    maxX = max(test_X)
    minX = min(test_X)
    maxY = max(test_Y)
    minY = min(test_Y)
    maxXY = max(maxX, maxY)

    parameterBounds = []
    parameterBounds.append([-maxXY, maxXY]) # seach bounds for a
    parameterBounds.append([-maxXY, maxXY]) # seach bounds for b
    parameterBounds.append([-maxXY, maxXY]) # seach bounds for c

    # "seed" the numpy random number generator for repeatable results
    result = differential_evolution(sumOfSquaredError, parameterBounds, seed=3)
    return result.x

# generate initial parameter values
geneticParameters = generate_Initial_Parameters()

# curve fit the test data
fittedParameters, pcov = curve_fit(func_power_law, test_X, test_Y, geneticParameters)

print('Parameters', fittedParameters)

modelPredictions = func_power_law(test_X, *fittedParameters) 

absError = modelPredictions - test_Y

SE = numpy.square(absError) # squared errors
MSE = numpy.mean(SE) # mean squared errors
RMSE = numpy.sqrt(MSE) # Root Mean Squared Error, RMSE
Rsquared = 1.0 - (numpy.var(absError) / numpy.var(test_Y))
print('RMSE:', RMSE)
print('R-squared:', Rsquared)

print()


##########################################################
# graphics output section
def ModelAndScatterPlot(graphWidth, graphHeight):
    f = plt.figure(figsize=(graphWidth/100.0, graphHeight/100.0), dpi=100)
    axes = f.add_subplot(111)

    # first the raw data as a scatter plot
    axes.plot(test_X, test_Y,  'D')

    # create data for the fitted equation plot
    xModel = numpy.linspace(min(test_X), max(test_X))
    yModel = func_power_law(xModel, *fittedParameters)

    # now the model as a line plot
    axes.plot(xModel, yModel)

    axes.set_xlabel('X Data') # X axis data label
    axes.set_ylabel('Y Data') # Y axis data label

    plt.show()
    plt.close('all') # clean up after using pyplot

graphWidth = 800
graphHeight = 600
ModelAndScatterPlot(graphWidth, graphHeight)