Scipy Optimize CurveFit вычисляет неправильные значения

Question

Мне интересно знать фазовый сдвиг между двумя типами синусоидальных волн. Для этого я пытаюсь подогнать каждую волну с помощью scipy.cuve_fit. Я подписался на этот пост . Однако я получаю отрицательные амплитуды, и фазовый сдвиг иногда выглядит как перенаправленные пи радианы.

Код, который я использую, приведен ниже:

def fit_sin_LD(t_LD, y_LD):

'''Fit sin to the input time sequence, and return fitting parameters "amp", "omega", "phase", "offset", "freq", "period" and "fitfunc"'''

ff = np.fft.fftfreq(len(t_LD), (t_LD[1]-t_LD[0]))   # assume uniform spacing
Fyy = abs(np.fft.fft(y_LD))
guess_freq = abs(ff[np.argmax(Fyy[1:])+1])   # excluding the zero frequency "peak", which is related to offset
guess_amp = np.std(y_LD) * 2.**0.5
guess_offset = np.mean(y_LD)
guess = np.array([guess_amp, 2.*np.pi*guess_freq, 0., guess_offset])

def sinfunc_LD(t_LD, A, w, p, c):
    return A * np.sin(w*t_LD + p) + c
#boundary=([0,-np.inf,-np.pi, 1.5],[0.8, +np.inf, np.pi, 2.5])
popt, pcov = scipy.optimize.curve_fit(sinfunc_LD, t_LD, y_LD, p0=guess, maxfev=3000) # with maxfev= number I can increase the number of iterations
A, w, p, c = popt
f          = w/(2.*np.pi)
fitfunc_LD = lambda t_LD: A*np.sin(w*t_LD + p) + c
fitted_LD  = fitfunc_LD(t_LD)
dic_LD = {"amp_LD": A, "omega_LD": w, "phase_LD": p, "offset_LD": c, "freq_LD": f, "period_LD": 1./f, "fitfunc_LD": fitted_LD, "maxcov_LD": np.max(pcov), "rawres_LD": (guess, popt, pcov)}
return dic_LD


def fit_sin_APD(t_APD, y_APD):

''' Fit sin to the input time sequence, and return fitting parameters "amp", "omega", "phase", "offset", "freq", "period" and "fitfunc" '''

ff = np.fft.fftfreq(len(t_APD), (t_APD[1]-t_APD[0]))   # assume uniform spacing
Fyy = abs(np.fft.fft(y_APD))
guess_freq = abs(ff[np.argmax(Fyy[1:])+1])   # excluding the zero frequency "peak", which is related to offset
guess_amp = np.std(y_APD) * 2.**0.5
guess_offset = np.mean(y_APD)
guess = np.array([guess_amp, 2.*np.pi*guess_freq, 0., guess_offset])

def sinfunc_APD(t_APD, A, w, p, c):
    return A * np.sin(w*t_APD + p) + c
#boundary=([0,0,-np.pi, 0.0],[np.inf, np.inf, np.pi, 0.7])
popt, pcov  = scipy.optimize.curve_fit(sinfunc_APD, t_APD, y_APD, p0=guess, maxfev=5000) # with maxfev= number I can increase the number of iterations
A, w, p, c  = popt
f           = w/(2.*np.pi)
fitfunc_APD = lambda t_APD: A*np.sin(w*t_APD + p) + c
fitted_APD  = fitfunc_APD(t_APD)
dic_APD     = {"amp_APD": A, "omega_APD": w, "phase_APD": p, "offset_APD": c, "freq_APD": f, "period_APD": 1./f, "fitfunc_APD": fitted_APD, "maxcov_APD": np.max(pcov), "rawres_APD": (guess, popt, pcov)}       
return dic_APD

Я не понимаю, почему функция curve_fit возвращает отрицательная амплитуда (что с точки зрения физики не имеет смысла). Я пробовал также устанавливать граничные условия как ** kwargs * с:

bounds=([0.0, -np.inf,-np.pi, 0.0],[+np.inf, +np.inf,-np.pi, +np.inf])

, но это дает более странный результат.

Я добавил изображение, показывающее эту разницу:

Кто-нибудь как побороть эту проблему с фазами и амплитудами?

Заранее спасибо

Answer 1

Здесь есть несколько вопросов, которые я не понимаю:

1. Нет необходимости определять функцию подгонки внутри «функции подгонки»
2. Нет необходимости определите его дважды, если единственная разница заключается в наименовании словаря. (Хотя я не понимаю, почему это должно называться по-разному в первую очередь)
3. Можно было бы подгонять частоту непосредственно вместо омега
4. При предварительном расчете подгоночных значений, непосредственно используйте заданная функция подгонки

В целом, я не понимаю, почему вторая подгонка должна потерпеть неудачу, и, используя некоторые общие данные c, это не так. Учитывая тот факт, что в физике амплитуда может быть сложной, у меня нет проблем с отрицательными результатами. Тем не менее, я понимаю смысл в ОП. Конечно, алгоритм подбора не знает физики, и математически нет проблем с отрицательной амплитудой. Это просто дает дополнительный сдвиг фазы пи. Следовательно, можно легко форсировать положительные амплитуды, заботясь о необходимом фазовом сдвиге. Я представил это здесь как возможный аргумент ключевого слова. Более того, я сократил это до одной функции подгонки с возможным «переименованием» выходных словарных ключей в качестве аргумента ключевого слова.

import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit


def sinfunc( t, A, f, p, c ):
    return A * np.sin( 2.0 * np.pi * f * t + p) + c

def fit_sin(t_APD, y_APD, addName="", posamp=False):

    ''' Fit sin to the input time sequence, and return fitting parameters "amp", "omega", "phase", "offset", "freq", "period" and "fitfunc" '''

    ff = np.fft.fftfreq( len( t_APD ), t_APD[1] - t_APD[0] ) # assume uniform spacing
    Fyy = abs( np.fft.fft( y_APD ) )
    guess_freq = abs( ff[np.argmax( Fyy[1:] ) + 1] )   # excluding the zero frequency "peak", which is related to offset
    guess_amp = np.std( y_APD ) * 2.**0.5
    guess_offset = np.mean( y_APD )
    guess = np.array( [ guess_amp, guess_freq, 0., guess_offset ] )

    popt, pcov  = curve_fit(sinfunc, t_APD, y_APD, p0=guess, maxfev=500) # with maxfev= number I can increase the number of iterations
    if popt[0] < 0 and posamp:
        popt[0] = -popt[0]
        popt[2] += np.pi 
        popt[2] = popt[2] % ( 2 * np.pi )
    A, f, p, c  = popt 
    fitted_APD  = sinfunc( t_APD, *popt )
    dic_APD     = {
            "amp{}".format(addName): A, 
            "omega{}".format(addName): 2.0 * np.pi * f, 
            "phase{}".format(addName): p, 
            "offset{}".format(addName): c, 
            "freq{}".format(addName): f, 
            "period{}".format(addName): 1.0 / f, 
            "fitfunc{}".format(addName): fitted_APD, 
            "maxcov{}".format(addName): np.max( pcov ), 
            "rawres{}".format(addName): ( guess, popt, pcov ) }
    return dic_APD

tl = np.linspace(0,1e-6, 150 )
sl1 = np.fromiter( (sinfunc(t, .18, 4998735, 3.6, 2.0 ) + .01 *( 1 - 2 * np.random.random() ) for t in tl ), np.float )
sl2 = np.fromiter( (sinfunc(t, .06, 4998735, 2.1, 0.4 ) + .01 *( 1 - 2 * np.random.random() ) for t in tl ), np.float )

ld = fit_sin(tl, sl1, addName="_ld" )
print ld["amp_ld"]
ld = fit_sin(tl, sl1, addName="_ld", posamp=True )
print ld["amp_ld"]
apd = fit_sin(tl, sl2 )

fig = plt.figure("1")
ax = fig.add_subplot( 1, 1, 1 )

ax.plot( tl, sl1, color="r" )
ax.plot( tl, ld["fitfunc_ld"], color="k", ls="--" )
ax.plot( tl, sl2, color="#50FF80" )
ax.plot( tl, apd["fitfunc"], color="k", ls="--" )

ax.grid()
plt.show()

Это дает мне:

-0.180108427200549
0.180108427200549

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