Scipy curve_fit дает неправильный ответ

Question

У меня есть осциллирующие данные, как показано на рисунке ниже, и я хочу подогнать их по синусоиде. Однако мой результат неверен.

Функция, которую я хочу подогнать под эту кривую:

def radius (z,phi, a0, k0,):

    Z = z.reshape(z.shape[0],1)

    k = np.array([k0,])
    a = np.array([a0,])
    r0 = 110
    rs  = r0 + np.sum(a*np.sin(k*Z +phi), axis=1)
    return rs

правильное решение может выглядеть следующим образом:

r_fit = radius(z, phi=np.pi/.8, a0=10,k0=0.017)
plt.plot(z, r,  label='data')
plt.plot(z,  r_fit, label='fitted curve')
plt.legend()

Мой результат, однако, из подгонки кривой выглядит так:

from scipy.optimize import curve_fit
popt, pcov = curve_fit(radius, xdata=z, ydata=r)

r_fit = radius(z, *popt)
plt.plot(z, r,  label='data')
plt.plot(z,  r_fit, label='fitted curve')
plt.legend()

Мои данные также следующие:

r = np.array([100.09061214, 100.17932773, 100.45526772, 102.27891728,
       113.12440802, 119.30644014, 119.86570527, 119.75184665,
       117.12160143, 101.55081608, 100.07280857, 100.12880236,
       100.39251753, 103.05404178, 117.15257288, 119.74048706,
       119.86955437, 119.37452005, 112.83384329, 101.0507198 ,
       100.05521567])

z = np.array([-407.90074345, -360.38004677, -312.99221012, -266.36934609,
       -224.36240585, -188.55933945, -155.21242348, -122.02778866,
        -87.84335638,  -47.0274899 ,    0.        ,   47.54559191,
         94.97469981,  141.33801462,  181.59490575,  215.77219256,
        248.95956379,  282.28027286,  318.16440024,  360.7246922 ,
        407.940799  ])

, поскольку моя функция просто представляет ряд Фурье, я также попробовал scipy.fftpack.fft (r), но не смог воспроизвести сигнал, близкий к тому, который ярассчитали БПФ.

Answer 1

Здесь приведен графический установщик Python с уравнением синуса и вашими данными с использованием модуля генетического алгоритма scipy.optimize diffrential Evolution для определения начальных оценок параметров для нелинейного решателя curve_fit. Этот модуль 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

r = numpy.array([100.09061214, 100.17932773, 100.45526772, 102.27891728,
       113.12440802, 119.30644014, 119.86570527, 119.75184665,
       117.12160143, 101.55081608, 100.07280857, 100.12880236,
       100.39251753, 103.05404178, 117.15257288, 119.74048706,
       119.86955437, 119.37452005, 112.83384329, 101.0507198 ,
       100.05521567])

z = numpy.array([-407.90074345, -360.38004677, -312.99221012, -266.36934609,
       -224.36240585, -188.55933945, -155.21242348, -122.02778866,
        -87.84335638,  -47.0274899 ,    0.        ,   47.54559191,
         94.97469981,  141.33801462,  181.59490575,  215.77219256,
        248.95956379,  282.28027286,  318.16440024,  360.7246922 ,
        407.940799  ])

 # rename data to match previous example code
xData = z
yData = r


def func (x, amplitude, center, width, offset): # equation sine[radians] + offset from zunzun.com
 return amplitude * numpy.sin(numpy.pi * (x - center) / width) + offset


# 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)

    diffY = maxY - minY
    diffX = maxX - minX

    parameterBounds = []
    parameterBounds.append([0.0, diffY]) # search bounds for amplitude
    parameterBounds.append([minX, maxX]) # search bounds for center
    parameterBounds.append([0.0, diffX]) # search bounds for width
    parameterBounds.append([minY, maxY]) # search bounds for offset

    # "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)

Answer 2

Проблема в том, что без предоставления первоначального предположения решение не может сходиться. Попробуйте добавить разумное начальное предположение:

p0 = [np.pi/.8, 10, 0.017]
popt, pcov = curve_fit(radius, xdata=z, ydata=r, p0=p0)

Обратите внимание, что если бы вы использовали один из других методов, таких как trf или dogbox, то без начального предположения это с большей вероятностью вернетОшибка RunTime из-за невозможности сближения параметров.