MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn moreOpportunities for recent engineering grads.

Apply Today
Asked by Kenneth on 8 Feb 2011

Hello all,

I need to fit a nonlinear model to several data sets simultaneously. The model has the same functional form for all sets, and the values of some model parameters are the same for all sets, but the value of at least one parameter is different for each set. For example, three exponential decay curves might have the same decay constant but a different amplitude for each data set. The global fit to all three curves would produce one decay constant and three amplitudes.

What is the best way to do this in MatLab? (Apologizing in advance for a novice question!)

ELELAB

*No products are associated with this question.*

Answer by Tom Lane on 11 Dec 2012

Accepted answer

Take a look at this example of a function that partitions the parameters into a set sharing values across all datasets and a set taking different values for different datasets. You may be able to modify this for your purposes.

function sharedparams t = (0:10)'; T = [t; t; t]; Y = 3 + [exp(-t/2); 2*exp(-t/2); 3*exp(-t/2)] + randn(33,1)/10; dsid = [ones(11,1); 2*ones(11,1); 3*ones(11,1)]; gscatter(T,Y,dsid)

X = [T dsid]; b = nlinfit(X,Y,@subfun,ones(1,5)) line(t,b(1)+b(2)*exp(-t/b(5)),'color','r') line(t,b(1)+b(3)*exp(-t/b(5)),'color','g') line(t,b(1)+b(4)*exp(-t/b(5)),'color','b')

function yfit = subfun(params,X) T = X(:,1); % time dsid = X(:,2); % dataset id A0 = params(1); % same A0 for all datasets A1 = params(2:4)'; % different A1 for each dataset tau = params(5); % same tau yfit = A0 + A1(dsid).*exp(-T/tau);

This uses nlinfit from the Statistics Toolbox, but you could use other fitting routines instead.

Tom Lane on 22 Jul 2013

I don't have your data, so I cannot try this, but I thought the output t from ode45 is supposed to match the input tspan. If that's the case, you can control its length.

Answer by j dr on 8 Feb 2011

Normaly you could identify both parameters in a fit using lsqcurvefit (least square fit) and a function of the type:

function F=myfun(x,xdata) F=x(1)*exp(-1*x(2)*xdata)

then

[x,resnorm] = lsqcurvefit(@myfun,StartingValues,xdata,ydata,LowerBound,UpperBound,options);

x will contain amplitude and decay rate that best match your data starting at StartingValue between LowerBound and UpperBound.

You can then appreciate that your 3 datasets have roughly the same decay rate and different amplitudes, OR:

If you know in advance which parameter will be shared, You can use lsqcurvefit with a function that uses a global variable as your shared parameter.

function F=myfunFixedAmp(x,xdata) global A F=A*exp(-1*x(1)*xdata)

or a different finction for a fixed decay rate

function F=myfunFixedDecay(x,xdata) global Alpha F=x(1)*exp(-1*Alpha*xdata)

Kenneth on 9 Feb 2011

Thanks, but I don't think this will do what I need. In the first example code you give (below) the amplitude is fixed and the decay constant is variable, but I need the amplitude to be a variable fitting parameter (which has the same value for all data sets) and the decay constant also to be a variable fitting parameter but it has a different value for each data set.

function F = myfunFixedAmp(x,xdata)

global A

F=A*exp(-1*x(1)*xdata)

Thanks,

ELELAB

Answer by Kenneth on 9 Feb 2011

(Wasn't sure if this should be a "comment" or "answer" so posting it as both)

Thanks, but I don't think this will do what I need. In the first example code you give (below) the amplitude is fixed and the decay constant is variable, but I need the amplitude to be a variable fitting parameter (which has the same value for all data sets) and the decay constant also to be a variable fitting parameter but it has a different value for each data set.

function F = myfunFixedAmp(x,xdata) global A F=A*exp(-1*x(1)*xdata)

So if there were three data sets to be fit, there would be four variable parameters: A (common to all three sets), and three decay constants, one for each set. One then minimizes the sum of the squared residuals of all three data sets.

Thanks again, ELELAB

Answer by j dr on 6 Apr 2011

Like I said in the top half of my post, if you use lsqcurvefit on each of the 3 curves independently, you would get:

A1,Alpha1 A2,Alpha2 A3,Alpha3

You could then verify if all "A"s are similar or all "Alpha"s are similar (by checking the standard deviation or the standard deviation divided by the mean to have a relative value of which of "A"s of "Alpha"s are more similar).

If this does not solve your problem I think you might need to restate it.

Answer by Richard Willey on 6 Apr 2011

Hi Kenneth

Coincidentially, I did a webinar a couple weeks back that uses lsqcurvefit to solve just this type of problem.

You can download the code from http://www.mathworks.com/matlabcentral/fileexchange/30869-fitting-with-matlab-statistics-optimization-and-curve-fitting

(You want the analysis.m example)

The webinar was titled "Fitting with MATLAB: Statistics, Optimization, and Curve Fitting". The recorded webinar should be up on the Statistics Toolbox product page shortly

Answer by DiRa on 11 Dec 2012

Dear all, it's been some time now since the last answer was posted, but I'm having almost the same problem now. I thought the last answer from Richard would help me, but it actually didn't. So here's my problem:

I have about 200 2D-datasets, where the x-scale shows the time(t)-dependence.

I need to fit a nonlinear exponetial decay function onto each one of these 200 datasets, that looks in principal like this:

y = y0 + A1*exp(-t/tau1) + A2*exp(-t/tau2) + A3*exp(-t/tau3)

The important thing is now that tau1,tau2 and tau3 must be the same for all 200 datasets and y0, A1, A2 and A3 may vary, but all of them have to be fitted by a lsqcurvefit at once. So the variables tau have to be constant for all datasets in the end, but is an unknown variable during the lsqcurvefit and should be optimized by the fit.

To specify it again: I want to do a lsqcurvefit for the variables y0,A1,A2,A3,tau1,tau2,tau3 for every single 2D dataset and in addition i need to do a global fit for tau1,tau2,tau3.

Is there a nice and fast way to calculate this kind of fit in matlab?

When I tried to fully understand the analysis.m i had the problem that it is not running after downloading all files because the is missing a needed function or script.

Thanks alot!

Answer by Seth DeLand on 11 Dec 2012

So if I understand correctly, you want to fit tau1, tau2, tau3 to all the data, but the other parameters are fit on an individual dataset basis? If so, I think this solution of breaking this up into a linear lsq problem embedded in a nonlinear lsq problem would work:

Identify tau1, tau2, tau3 for all of the data sets with lsqcurvefit. The other parameters y0, A1, A2, A3 can be identified on an individual basis inside the lsqcurvefit objective function (and this will be pretty fast because it's a linear least squares problem to solve for these).

Lsqcurvefit will be modifying the values for tau1, tau2, tau3, and then calling your objective function. Your objective function will:

1) Loop over the 200 datasets:

- solve a linear least-squares problem to identify y0,A1,A2,A3 for each dataset (this can be done with "x = A\b", which returns the least-squares soln for underdetermined systems: http://www.mathworks.com/help/matlab/ref/mldivide.html).
- calculate the value of your function at each point for each dataset: A*x

2) Combine all of the ydata from the 200 runs to be passed back to lsqcurvefit.

This approach may be kind of time consuming (it has to do 200 matrix solves each time the objective function is calculated), depending on the size of your datasets. But, it should be quicker and give you a better answer than trying to solve this as one large nonlinear curve fitting problem.

Answer by DiRa on 12 Dec 2012

Dear Seth and Tom,

I want to thank you both for your fast response. I used Toms suggestion since my script was already looking like that in some parts.

Since I needed some lower and upper bounds I then used the lsqcurvefit. After fixing some memory problems it then worked perfectly.

Thanks alot again!

Answer by Yi on 18 Apr 2013

@Kenneth: I understand precisely what your problem is and I am dealing with very similar situation: fitting 2D Gaussian function to hundreds of small pieces of ROI cut out from a large image, with the same global sigma, but different intensities and centroids. In fact in Igor, there is a built-in function called "global fit" that does the job. But honestly I don't know any elegant way to do the same in Matlab.

@Tom: Your solution will work for small set of data, like Kenneth's three exponential functions. But In my problem, I am dealing with greater than 500 sets of data at the same time, which would mean thousands of parameters to optimize at a time. But most of the parameters from different ROI don't "crosstalk", ie the Jacobian matrix is extremely spars. I am not sure if Matlab is smart enough to figure that out.

The workaround solution I can think of for large number of data sets like for my case is the E-M approach, which is to deal the common global parameter and other local parameters separately. But it can be slow too...

Answer by Julie on 14 May 2014

Dear fellow Matlab users I have a similar problem but can't seem to solve the issue with the above suggestions:

I am trying to use nlinfit to fit 14 sets of data to obtain a common exchange rate while a second scaling variable is also allowed to vary during the fit and should be different for each data set. I load a data file (f1) where the first row contains a 0 followed by my independent variables (24 of them) and each subsequent row has in column 1 a variable corresponding to the data in that row (ligand concentration) followed by the dependant variable (intensity) for the corresponding independent variables in row 1. I am getting the following error:

Error using nlinfit (line 122) Requires a vector second input argument.

Error in fit_lineshape_exchange (line 42) [beta,R]=nlinfit(x,y,@lineshape_exchange, beta0);

my function is below:

function f = lineshape_exchange(a)

%a(1) is k (the rate), a(2) is c (the weighting function allowing each

%spectrum to vary in intensity)

global K L v1 v2 Tp1 Tp2 w n

K=0.0000033;

Tp1=4.96;

Tp2=7.25;

v1=12.212;

v2=13.114;

f=0;

for j=1:n

t=1/(a(1*L*(1+K)));

p1=K*a(1)*L*t;

p2=1-p1;

W=v1-v2;

P=t*((1/Tp1*Tp2)-4*pi^2*w^2+pi^2*W^2)+(p1/Tp1)+(p2/Tp2);

Q=t*(2*pi*w-pi*W*(p1-p2));

R=2*pi*w*(1+t*((1/Tp1)+(1/Tp2)))+pi*W*t*((1/Tp1)-(1/Tp2)+pi*W*(p1-p2));

f=a(2)*(P*(1+t*(p2/Tp1 + p1/Tp2))+Q*R)/(P^2+R^2);

end

## 0 Comments