Как установить начальное значение для Curve_fit, чтобы найти лучшую оптимизацию, а не только локальную оптимизацию?

Question

Я пытаюсь подобрать степенную функцию и найти наиболее подходящий параметр.Тем не менее, я считаю, что если исходное предположение параметра отличается, вывод «наилучшего соответствия» будет другим.Если я не найду правильное первоначальное предположение, я смогу получить лучшую оптимизацию вместо локальной оптимизации.Есть ли способ найти ** подходящее начальное предположение ** ????.Мой код указан ниже.Пожалуйста, не стесняйтесь вносить любой вклад.Спасибо!

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)

, я смогу получить наилучшее соответствие

--------------------- Обновлено 10.7.2018 -------------------------------------------------------------------------------------------------------------------------

Поскольку мне нужно запустить тысячи или даже миллионов функции степенного закона, использование метода @James Phillips слишком дорого.Итак, какой метод подходит, кроме Curve_fit?такие как sklearn, np.linalg.lstsq и т. д.

Answer 1

Простого ответа нет: если бы он был, он был бы реализован в 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-fid как массив целых чисел , а затем выдается ошибка с отрицательными показателями.(И я полагаю, что цель 1.0 состояла в том, чтобы сделать его массивом с плавающей точкой? Для этого следует использовать параметр dtype=np.float.)

Answer 2

В дополнение к очень хорошим ответам из «Добро пожаловать в переполнение стека», что «не существует простого, универсального подхода и Джеймса Филлипса, что« дифференциальная эволюция часто помогает найти хорошие отправные точки (или даже хорошие решения!), Если они несколько медленнее, чем * 1001 »* ", позвольте мне дать отдельный ответ, который вам может пригодиться.

Во-первых, тот факт, что curve_fit() по умолчанию принимает любые значения параметров, является душераздирающей плохой идеей. Нет никакого возможного оправдания для такого поведения, и вы, и все остальные должны рассматривать тот факт, что для параметров существуют значения по умолчанию, как серьезную ошибку при реализации curve_fit() и делать вид, что этой ошибки не существует. НИКОГДА не считаю эти значения по умолчанию разумными.

Из простого графика данных должно быть очевидно, что a=1, b=1, c=1 являются очень, очень плохими начальными значениями. Функция затухает, поэтому b < 0. На самом деле, если бы вы начали с a=1, b=-1, c=1, вы бынашли правильное решение.

Возможно, это также помогло быnd по параметрам.Даже установление границ c (-100, 100), возможно, помогло.Как и в случае знака b, я думаю, что вы могли видеть эту границу из простого графика данных.Когда я пытаюсь сделать это для вашей проблемы, границы для c не помогают, если начальное значение равно b=1, но для b=0 или b=-5.

Что еще более важно, хотя вы печатаетеНа графике лучше всего подходят параметры popt. Вы не печатаете неопределенности или корреляции между переменными, содержащимися в pcov, и, таким образом, ваша интерпретация результатов является неполной.Если бы вы посмотрели на эти значения, вы бы увидели, что начиная с b=1 приводит не только к плохим значениям, но также к огромной неопределенности параметров и очень, очень высокой корреляции.Это припадок говорит вам, что он нашел плохое решение.К сожалению, возврат pcov из curve_fit не очень легко распаковать.

Позвольте мне порекомендовать 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

, что несколько лучше.

Короче, начальные значения всегда имеют значение, иВ результате получаются не только наиболее подходящие значения, но и неопределенности и корреляции.

Answer 3

Вот пример кода с использованием генетического алгоритма scipy.optimize.differential_evolution с вашими данными и уравнением.Этот модуль scipy использует алгоритм Latin Hypercube для обеспечения тщательного поиска пространства параметров и поэтому требует границ, в которых искать - в этом примере эти границы основаны на максимальных и минимальных значениях данных.Для других проблем вам может потребоваться указать другие границы поиска, если вы знаете, какой диапазон значений параметров следует ожидать.

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)