如何设置 curve_fit 的初始值以找到最佳优化，而不仅仅是局部优化？
how to set up the initial value for curve_fit to find the best optimizing, not just local optimizing?
我正在尝试拟合幂律函数，以便找到最合适的参数。但是，我发现如果参数的初始猜测不同，"best fit" 输出就会不同。除非我找到正确的初始猜测，否则我可以获得最佳优化，而不是局部优化。有没有办法找到**合适的初始猜测**？？？？。下面列出了我的代码。请随时输入任何信息。谢谢！
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 年 7 月 10 日更新---------- ---------------------------------------------- ---------------------------------------------- ------------
由于我需要 运行 数千甚至 数百万 的幂律函数，使用@James Phillips 的方法太昂贵了。那么除了curve_fit还有什么方法合适呢？比如sklearn，np.linalg.lstsq等
这是使用 scipy.optimize.differential_evolution 遗传算法的示例代码，包含您的数据和方程。此 scipy 模块使用拉丁超立方体算法来确保对参数 space 进行彻底搜索，因此需要搜索范围 - 在本例中，这些范围基于数据的最大值和最小值。对于其他问题，如果您知道预期的参数值范围，则可能需要提供不同的搜索范围。
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)
没有简单的答案：如果有，它会在 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 参数的目的。）
除了来自 Welcome to Stack Overflow 的非常好的答案之外，"there is no easy, universal approach and James Phillips that "差异进化经常
如果比 curve_fit()" 慢一些，有助于找到好的起点（甚至好的解决方案！），请允许我给出一个单独的答案，您可能会觉得有帮助。
首先，curve_fit() 默认为任何参数值这一事实是令人心碎的坏主意。这种行为没有可能的理由，你和其他人应该把参数有默认值这一事实视为 curve_fit() 实现中的严重错误，并假装这个错误不存在。 从不相信这些默认值是合理的。
从一个简单的数据图中可以明显看出，a=1, b=1, c=1 是非常非常糟糕的初始值。函数衰减，所以b < 0。事实上，如果您从 a=1, b=-1, c=1 开始，您就会找到正确的解决方案。
它可能还有助于对参数设置合理的界限。即使设置 c of (-100, 100) 的边界也可能有所帮助。与 b 的符号一样，我认为您可以从简单的数据图中看到该边界。当我为你的问题尝试这个时，如果初始值为 b=1，c 上的界限没有帮助，但它对 b=0 或 b=-5.
有帮助
更重要的是，虽然您在图中打印了最合适的参数 popt，但您没有打印 pcov 中变量之间的不确定性或相关性，因此您对结果的解释不完整。如果您查看过这些值，您会发现以 b=1 开头不仅会导致错误的值，还会导致参数的巨大不确定性和非常非常高的相关性。这恰好告诉你它找到了一个糟糕的解决方案。不幸的是，来自 curve_fit 的 return 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
哪个好一点。
简而言之，初始值总是很重要，结果不仅是最佳拟合值，还包括不确定性和相关性。
我正在尝试拟合幂律函数，以便找到最合适的参数。但是，我发现如果参数的初始猜测不同，"best fit" 输出就会不同。除非我找到正确的初始猜测，否则我可以获得最佳优化，而不是局部优化。有没有办法找到**合适的初始猜测**？？？？。下面列出了我的代码。请随时输入任何信息。谢谢！
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 年 7 月 10 日更新---------- ---------------------------------------------- ---------------------------------------------- ------------
由于我需要 运行 数千甚至 数百万 的幂律函数，使用@James Phillips 的方法太昂贵了。那么除了curve_fit还有什么方法合适呢？比如sklearn，np.linalg.lstsq等
这是使用 scipy.optimize.differential_evolution 遗传算法的示例代码，包含您的数据和方程。此 scipy 模块使用拉丁超立方体算法来确保对参数 space 进行彻底搜索，因此需要搜索范围 - 在本例中，这些范围基于数据的最大值和最小值。对于其他问题，如果您知道预期的参数值范围，则可能需要提供不同的搜索范围。
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)
没有简单的答案：如果有，它会在 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 参数的目的。）
除了来自 Welcome to Stack Overflow 的非常好的答案之外，"there is no easy, universal approach and James Phillips that "差异进化经常
如果比 curve_fit()" 慢一些，有助于找到好的起点（甚至好的解决方案！），请允许我给出一个单独的答案，您可能会觉得有帮助。
首先，curve_fit() 默认为任何参数值这一事实是令人心碎的坏主意。这种行为没有可能的理由，你和其他人应该把参数有默认值这一事实视为 curve_fit() 实现中的严重错误，并假装这个错误不存在。 从不相信这些默认值是合理的。
从一个简单的数据图中可以明显看出，a=1, b=1, c=1 是非常非常糟糕的初始值。函数衰减，所以b < 0。事实上，如果您从 a=1, b=-1, c=1 开始，您就会找到正确的解决方案。
它可能还有助于对参数设置合理的界限。即使设置 c of (-100, 100) 的边界也可能有所帮助。与 b 的符号一样，我认为您可以从简单的数据图中看到该边界。当我为你的问题尝试这个时，如果初始值为 b=1，c 上的界限没有帮助，但它对 b=0 或 b=-5.
更重要的是，虽然您在图中打印了最合适的参数 popt，但您没有打印 pcov 中变量之间的不确定性或相关性，因此您对结果的解释不完整。如果您查看过这些值，您会发现以 b=1 开头不仅会导致错误的值，还会导致参数的巨大不确定性和非常非常高的相关性。这恰好告诉你它找到了一个糟糕的解决方案。不幸的是，来自 curve_fit 的 return 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
哪个好一点。
简而言之，初始值总是很重要，结果不仅是最佳拟合值，还包括不确定性和相关性。