在python中用固定参数拟合函数和

Question

我要拟合的信号是多个正弦函数（和噪声）的叠加，我想同时适用于所有频率。这里有一个示例数据文件，由两个频率240d^-1和261.8181d^-1生成： https://owncloud.gwdg.de/index.php/s/JZQTJ3VMYZH8qNB和plot of the time series (excerpt)

到目前为止，我可以拟合一个接一个的正弦函数，同时保持频率固定在一个值上。我从周期图中得到频率，最后我对拟合的振幅和相位感兴趣。在

import numpy as np
from scipy import optimize
import bottleneck as bn

def f_sinus0(x,a,b,c,d):
    return a*np.sin(b*x+c)+d

def fit_single(t, flux, flux_err, freq_model, c0 = 0.):

    # initial guess for the parameter
    d0 = bn.nanmean(flux)
    a0 = 3*np.std(flux)/np.sqrt(2.)

    # fit function with fixed frequency "freq_model"
    popt, pcov = optimize.curve_fit(lambda x, a, c, d:
        f_sinus0(x, a, freq_model*2*np.pi, c, d),
        t, flux, sigma = flux_err, p0 = (a0,c0,d0),
        bounds=([a0-0.5*abs(a0),-np.inf,d0-0.25*abs(d0)],
        [a0+0.5*abs(a0),np.inf,d0+0.25*abs(d0)]),
        absolute_sigma=True)
    perr = np.sqrt(np.diag(pcov))

    return popt, perr

filename = 'data-test.csv'

data = np.loadtxt(filename)
time = data[0]
flux = data[1]
flux_err = data[2]

freq_model = 260 #d^-1

popt, perr = fit_single(time, flux, flux_err, freq_model, c0 = 0.)

现在我想同时拟合两个频率。我定义了一个函数，它根据输入参数列表的长度返回拟合函数的和

^{pr2}$

进行拟合

def fit_multiple(t, flux, flux_err, guess):
    popt, pcov = optimize.curve_fit(
        f_multiple_sin, t, flux, sigma=flux_err, p0=guess,
        bounds=(guess-np.multiply(guess,0.1),guess+np.multiply(guess,0.1)),
        absolute_sigma=True
        )

    perr = np.sqrt(np.diag(pcov))

    return popt, perr

guess = [4.50148944e-03, 2.40000040e+02, 3.01766641e-03, 8.99996136e-01, 3.14546648e-03, 2.61818207e+02, 2.94282247e-03, 5.56770657e-06]
popt, perr = fit_multiple(time, flux, flux_err, guess)

使用单个拟合的结果作为初始参数guess = [amplitude1, frequency1, phase1, offset1, amplitude2,...]

但是我怎样才能拟合出多个正弦函数，每个函数都有一个固定的频率？在这种情况下，lambda方法对我来说似乎不是那么直接。在

Answer 1

这是一个使用scipy.optimize.leastsq的解决方案，它给了我更多的自由。不过，在错误评估时，你必须小心。另一方面，在参数数目方面，它没有curve_fit那么严格。 在这个解决方案中，我基本上符合三个列表，振幅、频率和相位。At似乎很方便将它按这样的顺序传递给函数。 最后你可以固定任何频率子集。不过，我的印象是，收敛对起始参数非常敏感。在

import matplotlib.pyplot as plt
import numpy as np
import scipy.optimize as so


def multisine(x, ampList, freqList, phaseList):
    assert len( ampList ) == len( freqList )
    assert len( ampList ) == len( phaseList )
    out=0
    for a, f, p in zip( ampList, freqList, phaseList ):
        out += a * np.sin( x * f + p )
    return out


### FixedFrequencies is a list of values and positions in the list to pass to multisine....remember counting from zero
def multisine_fixed_fs(x, params, n, FixedFrequencies=None):
    if FixedFrequencies is None:
        assert len( params ) == 3 *  n
        ampList = params[ : n]
        freqList = params[ n : 2* n] 
        phaseList = params[ 2 * n : ]
    else:
        assert len( params ) + len( FixedFrequencies ) == 3 *  n
        ampList = params[ : n]
        freqList = list(params[ n : -n ])
        phaseList = params[ -n : ]
        sortedList = sorted( list(FixedFrequencies), key=lambda x: x[-1] )
        for fixed in sortedList:
            freqList.insert(fixed[-1], fixed[0] )

    return multisine(x, ampList, freqList, phaseList)


def residuals(params, data, n, FixedFrequencies=None):
    xList, yList = zip( *data )
    thyList = [ multisine_fixed_fs( x, params, n , FixedFrequencies=FixedFrequencies ) for x in xList ]
    d = [ y1- y2 for y1, y2 in zip( yList, thyList ) ]
    return d



xList = np.linspace( 0, 100, 100 )
yList = np.fromiter( ( multisine(x, [ 1, .3 ], [ .4, .42 ],[ 0, .1] ) for x in xList ), np.float )

data = zip( xList, yList )

fit, err = so.leastsq( residuals,  x0=[ 1.2, .32 ] + [ .42, .43 ] + [ 0.1, 0.12 ], args=( data, 2 ) )
print fit

fit, err = so.leastsq( residuals,  x0=[ 1.2, .32 ] + [ .42 ] + [ 0.1, 0.12 ], args=( data, 2 , [ [ .45, 1 ] ]) )
print fit
y2List = np.fromiter( ( multisine(x, [ fit[0], fit[1] ], [ fit[2], .45 ],[ fit[-2], fit[-1] ] ) for x in xList ), np.float )

fit, err = so.leastsq( residuals,  x0=[ 1.2, .32 ]  + [ 0.1, 0.12 ], args=( data, 2 , [ [ .39, 0 ],[ .45, 1 ] ]) )
print fit
y3List = np.fromiter( ( multisine(x, [ fit[0], fit[1] ], [ .39, .45 ],[ fit[-2], fit[-1] ] ) for x in xList ), np.float )

fig = plt.figure(1)
ax = fig.add_subplot( 1, 1, 1 )
ax.plot(xList,yList)
ax.plot(xList,y2List)
ax.plot(xList,y3List)

plt.show()

提供：

^{pr2}$