我正在尝试拟合幂律函数，以便找到最合适的参数。但是，我发现如果参数的初始猜测不同，则“最适合”输出是不同的。除非我找到正确的初始猜测，否则我可以获得最佳优化，而不是局部优化。有没有办法找到**合适的初始猜测**？？？？。下面列出了我的代码。请随时输入任何信息。谢谢!

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 模块使用拉丁超立方体算法来确保对参数空间进行彻底搜索，因此需要在搜索范围内进行搜索——在本例中，这些范围基于数据的最大值和最小值。对于其他问题，如果您知道预期的参数值范围，则可能需要提供不同的搜索范围。

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)