2

I am trying to write a code that would fit multiple (>100) data sets with the same parameters.

For a simple example:

Suppose I have 4 different data sets that are supposed to approximately follow the function y(x)=A*Sin(bx+c)+dx. They all have the same parameters a,b,and d. The only parameter that is different is c (phase), but for each data set, I know exactly what c should be.

Now I want to fit all 4 of those functions simultaneously and get one set of fit parameters a,b and d out.

Here I wrote a code that fits each of the functions separately and gives me 4 different sets of parameters, which is not what I am looking for.

Thank you in advance for any help.

from scipy import optimize
import matplotlib
from matplotlib import pyplot as plt
from scipy import *

Mat_ydata=[]
Mat_angle=[]
Mat_xdata=[]

defining a function for fitting

def sin_func(x,a,b,c,d):
   return a*sin(b*x+c)+d*x

Generating some data: 4 sine functions with the same amplitude,frequency and linear offset, but each has different phase, that we know already. Also each is defined over different range. The x and y corrdinates of the functions are stored in the lists of lists Mat_ydata and Mat_xdata. The phase is stored in the list Mat_phase

for phase in range(0,181,60):
    angle=float(phase)/180*pi
    num_points=phase+50
    xdata=linspace(num_points/10,num_points/5, num_points)
    ydata=2.11*sin(pi/2*xdata+angle)+3.16*xdata+2*((0.5-rand(num_points))*exp(2*rand(num_points)**2))
    Mat_angle.append(angle)
    Mat_ydata.append(ydata)
    Mat_xdata.append(xdata)

Here is fitting. I am fitting each function separately right now and therefore getting different values for phase, amplitude and linear offset for each of them. I would like to fit them all together. In the end I also plot data and fits.

for i in range(0,len(Mat_angle)):
    fitfunc = lambda p, x: sin_func(x,p[0],p[1],Mat_angle[i],p[2])
    errfunc = lambda p, x, y: fitfunc(p, x) - y 
    p0 = [2,1.5,3] 
    p1, success = optimize.leastsq(errfunc, p0[:], args=(Mat_xdata[i], Mat_ydata[i]))
    print p1
    plt.plot(Mat_xdata[i],Mat_ydata[i],'o',Mat_xdata[i],fitfunc(p1,Mat_xdata[i]),'-')
plt.show()
PumpkinPie
  • 111
  • 1
  • 6

1 Answers1

4

I have actually figured it out in the last couple of days. I'll provide the code in case it may be of interest to anyone. I also found that fitting sine functions is pretty hard, so I changed my fake data to Lorentzians.

Importing modules and generating fake data, that will be kept in the the lists of lists:

from scipy import optimize
import matplotlib
from matplotlib import pyplot as plt
from scipy import *
import numpy as np


Mat_ydata=[]
Mat_angle=[]
Mat_xdata=[]

for c in range(0,100,20):
    num_points=c+100
    xdata=linspace(0,num_points/2, num_points)
    ydata=5.1/((xdata-c)**2+2.1**2)+0.05*((0.5rand(num_points))*exp(2*rand(num_points)**2))
    Mat_angle.append(c)
    Mat_ydata.append(ydata)
    Mat_xdata.append(xdata)

Defining the fitting function and the error function:

def lor_func(x,c,par):
    a,b,d=par
    return a/((x-c)**2+b**2)

def err (p,c,x,y):
    return lor_func(x,c,p)-y

def err_global(p,Mat_a,Mat_x,Mat_y):
    err0=[]
    for i in range(0, len(Mat_a)):
        errc=err(p,Mat_a[i],Mat_x[i],Mat_y[i])
        err1=np.concatenate((err0,errc))
        err0=err1
    return err0

Actual fitting and displaying results:

p_global=[1,1,1]
p_best,success=optimize.leastsq(err_global, p_global,args=(Mat_angle,Mat_xdata,Mat_ydata),maxfev=40000)
toplot=[]
for i in range(0,len(Mat_angle)):
    toplot.append(lor_func(Mat_xdata[i],Mat_angle[i],p_best))
err_toplot=err_global(p_best,Mat_angle,Mat_xdata,Mat_ydata)
print p_best

for i in range(0,len(Mat_angle)):
    plt.plot(Mat_xdata[i],Mat_ydata[i],'o',Mat_xdata[i],toplot[i],'-')
plt.show()

And here is what it displays:

enter image description here

PumpkinPie
  • 111
  • 1
  • 6