Выбор границ Python curve_fit и начальное условие влияют на результат

У меня есть набор данных, описываемый двумя свободными параметрами, которые я хочу определить с помощью optimalization.curve_fit. Модель определяется следующим образом

def func(x, a, b,):
    return a*x*np.sqrt(1-b*x)

И подходящая часть как

popt, pcov = opt.curve_fit(f = func, xdata = x_data, ydata= y_data, p0 
= init_guess, bounds = ([a_min, b_min], [a_max, b_max]))

Результат решений для a и b довольно сильно зависит от моего выбора init_guess, т. е. начального предположения, а также от выбора границ. Есть ли способ решить это?


python model-fitting curve-fitting boundary
person Rob Schneider    schedule 20.06.2017    source источник
comment
То, что выбор границ влияет на решение, не так уж удивительно, равно как и первоначальные догадки. Я предполагаю, что лучший способ справиться с этим - это i) выбирать границы только тогда, когда вы действительно уверены в них (например, в вашем примере границы для b должны быть выбраны таким образом, что определяется sqrt) и ii) чтобы посмотрите на качество подгонки, построив график, рассчитав остатки и т. д. Добавление некоторых данных также может помочь, чтобы люди могли немного поиграть... :)   -  person Cleb    schedule 20.06.2017

Ответы (1)


arrow_upward
4
arrow_downward

Авторы модуля Python scipy включили генетический алгоритм Differential Evolution в код оптимизации scipy как модуль scipy.optimize.differential_evolution. Этот модуль можно использовать для стохастического поиска начальных значений параметров для нелинейной регрессии.

Вот пример кода из RamanSpectroscopyFit, который использует генетический алгоритм scipy для начальной оценки параметров для подбора данных рамановской спектроскопии:

import numpy as np
import pickle # for loading pickled test data
import matplotlib
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import warnings

from scipy.optimize import differential_evolution


# Double Lorentzian peak function
# bounds on parameters are set in generate_Initial_Parameters() below
def double_Lorentz(x, a, b, A, w, x_0, A1, w1, x_01):
    return a*x+b+(2*A/np.pi)*(w/(4*(x-x_0)**2 + w**2))+(2*A1/np.pi)*(w1/(4*(x-x_01)**2 + w1**2))


# function for genetic algorithm to minimize (sum of squared error)
# bounds on parameters are set in generate_Initial_Parameters() below
def sumOfSquaredError(parameterTuple):
    warnings.filterwarnings("ignore") # do not print warnings by genetic algorithm
    return np.sum((yData - double_Lorentz(xData, *parameterTuple)) ** 2)


def generate_Initial_Parameters():
    # min and max used for bounds
    maxX = max(xData)
    minX = min(xData)
    maxY = max(yData)
    minY = min(yData)

    parameterBounds = []
    parameterBounds.append([-1.0, 1.0]) # parameter bounds for a
    parameterBounds.append([maxY/-2.0, maxY/2.0]) # parameter bounds for b
    parameterBounds.append([0.0, maxY*100.0]) # parameter bounds for A
    parameterBounds.append([0.0, maxY/2.0]) # parameter bounds for w
    parameterBounds.append([minX, maxX]) # parameter bounds for x_0
    parameterBounds.append([0.0, maxY*100.0]) # parameter bounds for A1
    parameterBounds.append([0.0, maxY/2.0]) # parameter bounds for w1
    parameterBounds.append([minX, maxX]) # parameter bounds for x_01

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



# load the pickled test data from original Raman spectroscopy
data = pickle.load(open('data.pkl', 'rb'))
xData = data[0]
yData = data[1]

# generate initial parameter values
initialParameters = generate_Initial_Parameters()

# curve fit the test data
fittedParameters, niepewnosci = curve_fit(double_Lorentz, xData, yData, initialParameters)

# create values for display of fitted peak function
a, b, A, w, x_0, A1, w1, x_01 = fittedParameters
y_fit = double_Lorentz(xData, a, b, A, w, x_0, A1, w1, x_01)

plt.plot(xData, yData) # plot the raw data
plt.plot(xData, y_fit) # plot the equation using the fitted parameters
plt.show()

print(fittedParameters)
person James Phillips    schedule 20.06.2017