Fitting a sum of functions with fixed parameter in python

Question

The signal I want to fit is a superposition of multiple sine-functions (and noise) and I want to fit for all frequencies simultaneously. Here an example data file, generated with two frequencies of 240d^-1 and 261.8181d^-1: https://owncloud.gwdg.de/index.php/s/JZQTJ3VMYZH8qNB and plot of the time series (excerpt)

So far I can fit one sine-function after the other, while keeping the frequency fixed to a value. I get the frequency from e.g. a periodogram and in the end I am interested in amplitude and phase of the fit.

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

Now I want to fit both frequencies simultaneously. I defined a function that returns a sum of fitting-functions, depending on the length of the input-parameter-list like this

def f_multiple_sin(x, *params):
    y = np.zeros_like(x)
    for i in range(0, len(params), 4): #4=amplitude, freq, phase, offset
        amplitude = params[i]
        freq = params[i+1]
        phase = params[i+2]
        offset = params[i+3]
        y = y + amplitude*np.sin(np.multiply(freq, x)+phase)+offset
    return y

Performing the fit

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)

using the results from the individual fits as initial parameters guess = [amplitude1, frequency1, phase1, offset1, amplitude2,...]

But how can I fit multiple sine-functions, each with a fixed frequency? The lambda approach seems not so straight forward to me in this case.

Answer 1

This is a solution using scipy.optimize.leastsq which gives me more freedom. On error evaluation you have to take some care, though. On the other hand it is not as strict as curve_fit concerning the number of parameters. In this solution I fit basically three lists, the amplitudes, frequencies, and phases. At seemed convenient to pass it sorted like this to the function. At the end you can fix any subset of frequencies. I had the impression, though, that convergence is very sensitive to starting parameters.

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

Providing:

>> [ 1.00000006e+00   2.99999889e-01   3.99999999e-01   4.20000009e-01 1.47117910e-07   6.38318486e+00 ]
>> [ 1.12714624  0.12278804  0.40198029  0.08039605 -1.08564396 ]
>> [ 1.05124097 -0.32600116  0.6633511   1.18400026 ]