Почему Python's curve_fit не заканчивает оптимизацию? - PullRequest
Новая
       36

Почему Python's curve_fit не заканчивает оптимизацию?

1 голос
/ 10 июня 2019

Мне нужно найти два параметра уравнения, которые наилучшим образом соответствуют заданным значениям x и y.

Я использую Python 3 с Numpy и Scipy. 

from scipy.optimize import curve_fit

def func(dx, d50, p):
    return (1 / (1 + ((d50 / dx) ** p)))

xdata = [280, 150, 75, 45, 38, 20, 10, 5.1, 2.6]
ydata = [99.57592773, 95.53773499, 81.14313507, 67.08183289, 62.93716431, 49.961483, 37.80876923, 24.53152657, 13.2219696]

# curve fit:
popt, pcov = curve_fit(func, xdata, ydata)
print(popt)

I expect a d50 ~ 20 and a p > 0.

Но Python отправляет мне: 

[0.00221498 1.60291553]

> /usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:4:
> RuntimeWarning: invalid value encountered in power

после удаления cwd из sys.path.

1 Ответ

0 голосов
/ 10 июня 2019

Мне не удалось получить хорошее соответствие вашим данным, используя уравнение в вашем посте. Мой поиск по уравнению обнаружил, что стандартное уравнение пика Вейбулла, «a * exp (-0,5 * pow (log (x / b) / c, 2.0))», дает RMSE = 1,619 и R-квадрат = 0,997 для параметров a = 103.1533969 , b = 498,93546398 и с = 2,67321918, как показано ниже. Я включил графического установщика Python, использующего это уравнение, и стандартный модуль генетического алгоритма scipy diff_evolution, чтобы найти начальные оценки параметров для curve_fit (), этот модуль scipy использует алгоритм Latin Hypercube для обеспечения тщательного поиска пространства параметров, и этот алгоритм требует границ в пределах который искать. В этом примере границы поиска выводятся из данных. Гораздо проще определить диапазоны для начальных оценок параметров, чем найти конкретные значения.

plot

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


xData = [280, 150, 75, 45, 38, 20, 10, 5.1, 2.6]
yData = [99.57592773, 95.53773499, 81.14313507, 67.08183289, 62.93716431, 49.961483, 37.80876923, 24.53152657, 13.2219696]


def func(x, a, b, c): # Peak_WeibullPeak_model from zunzun.com
    return a * numpy.exp(-0.5 * numpy.power(numpy.log(x/b) / c, 2.0))


# 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(xData, *parameterTuple)
    return numpy.sum((yData - val) ** 2.0)


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

    minData = min(minX, minY)
    maxData = max(maxY, maxX)

    parameterBounds = []
    parameterBounds.append([minData, maxData]) # search bounds for a
    parameterBounds.append([minData, maxData]) # search bounds for b
    parameterBounds.append([minData, maxData]) # search bounds for c

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

# by default, differential_evolution completes by calling curve_fit() using parameter bounds
geneticParameters = generate_Initial_Parameters()

# now call curve_fit without passing bounds from the genetic algorithm,
# just in case the best fit parameters are aoutside those bounds
fittedParameters, pcov = curve_fit(func, xData, yData, geneticParameters)
print('Fitted parameters:', fittedParameters)
print()

modelPredictions = func(xData, *fittedParameters) 

absError = modelPredictions - yData

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(yData))

print()
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(xData, yData,  'D')

    # create data for the fitted equation plot
    xModel = numpy.linspace(min(xData), max(xData))
    yModel = func(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)
...