Simple nonlinear least squares curve fitting in MATLAB

December 6th, 2013 | Categories: math software, matlab | Tags:

A question I get asked a lot is ‘How can I do nonlinear least squares curve fitting in X?’ where X might be MATLAB, Mathematica or a whole host of alternatives.  Since this is such a common query, I thought I’d write up how to do it for a very simple problem in several systems that I’m interested in

This is the MATLAB version. For other versions,see the list below

The problem

xdata = -2,-1.64,-1.33,-0.7,0,0.45,1.2,1.64,2.32,2.9
ydata = 0.699369,0.700462,0.695354,1.03905,1.97389,2.41143,1.91091,0.919576,-0.730975,-1.42001

and you’d like to fit the function using nonlinear least squares.  You’re starting guesses for the parameters are p1=1 and P2=0.2

For now, we are primarily interested in the following results:

• The fit parameters
• Sum of squared residuals

Future updates of these posts will show how to get other results such as confidence intervals. Let me know what you are most interested in.

How you proceed depends on which toolboxes you have. Contrary to popular belief, you don’t need the Curve Fitting toolbox to do curve fitting…particularly when the fit in question is as basic as this. Out of the 90+ toolboxes sold by The Mathworks, I’ve only been able to look through the subset I have access to so I may have missed some alternative solutions.

Pure MATLAB solution (No toolboxes)

In order to perform nonlinear least squares curve fitting, you need to minimise the squares of the residuals. This means you need a minimisation routine.  Basic MATLAB comes with the fminsearch function which is based on the Nelder-Mead simplex method.  For this particular problem, it works OK but will not be suitable for more complex fitting problems.  Here’s the code

format compact
format long
xdata = [-2,-1.64,-1.33,-0.7,0,0.45,1.2,1.64,2.32,2.9];
ydata = [0.699369,0.700462,0.695354,1.03905,1.97389,2.41143,1.91091,0.919576,-0.730975,-1.42001];

%Function to calculate the sum of residuals for a given p1 and p2
fun = @(p) sum((ydata - (p(1)*cos(p(2)*xdata)+p(2)*sin(p(1)*xdata))).^2);

%starting guess
pguess = [1,0.2];

%optimise
[p,fminres] = fminsearch(fun,pguess)

This gives the following results

p =
1.881831115804464 0.700242006994123
fminres =
0.053812720914713

All we get here are the parameters and the sum of squares of the residuals.  If you want more information such as 95% confidence intervals, you’ll have a lot more hand-coding to do. Although fminsearch works fine in this instance, it soon runs out of steam for more complex problems.

MATLAB with optimisation toolbox

With respect to this problem, the optimisation toolbox gives you two main advantages over pure MATLAB.  The first is that better optimisation routines are available so more complex problems (such as those with constraints) can be solved and in less time.  The second is the provision of the lsqcurvefit function which is specifically designed to solve curve fitting problems.  Here’s the code

format long
format compact
xdata = [-2,-1.64,-1.33,-0.7,0,0.45,1.2,1.64,2.32,2.9];
ydata = [0.699369,0.700462,0.695354,1.03905,1.97389,2.41143,1.91091,0.919576,-0.730975,-1.42001];

%Note that we don't have to explicitly calculate residuals
fun = @(p,xdata) p(1)*cos(p(2)*xdata)+p(2)*sin(p(1)*xdata);

%starting guess
pguess = [1,0.2];

[p,fminres] = lsqcurvefit(fun,pguess,xdata,ydata)

This gives the results

p =
1.881848414551983   0.700229137656802
fminres =
0.053812696487326

MATLAB with statistics toolbox

There are two interfaces I know of in the stats toolbox and both of them give a lot of information about the fit. The problem set up is the same in both cases

%set up for both fit commands in the stats toolbox
xdata = [-2,-1.64,-1.33,-0.7,0,0.45,1.2,1.64,2.32,2.9];
ydata = [0.699369,0.700462,0.695354,1.03905,1.97389,2.41143,1.91091,0.919576,-0.730975,-1.42001];

fun = @(p,xdata) p(1)*cos(p(2)*xdata)+p(2)*sin(p(1)*xdata);

pguess = [1,0.2];

Method 1 makes use of NonLinearModel.fit

mdl = NonLinearModel.fit(xdata,ydata,fun,pguess)

The returned object is a NonLinearModel object

class(mdl)
ans =
NonLinearModel

which contains all sorts of useful stuff

mdl =
Nonlinear regression model:
y ~ p1*cos(p2*xdata) + p2*sin(p1*xdata)

Estimated Coefficients:
Estimate             SE
p1      1.8818508110535      0.027430139389359
p2    0.700229815076442    0.00915260662357553

tStat               pValue
p1    68.6052223191956    2.26832562501304e-12
p2    76.5060538352836    9.49546284187105e-13

Number of observations: 10, Error degrees of freedom: 8
Root Mean Squared Error: 0.082
R-Squared: 0.996,  Adjusted R-Squared 0.995
F-statistic vs. zero model: 1.43e+03, p-value = 6.04e-11

If you don’t need such heavyweight infrastructure, you can make use of the statistic toolbox’s nlinfit function

[p,R,J,CovB,MSE,ErrorModelInfo] = nlinfit(xdata,ydata,fun,pguess);

Along with our parameters (p) this also provides the residuals (R), Jacobian (J), Estimated variance-covariance matrix (CovB), Mean Squared Error (MSE) and a structure containing information about the error model fit (ErrorModelInfo).

Both nlinfit and NonLinearModel.fit use the same minimisation algorithm which is based on Levenberg-Marquardt

MATLAB with Symbolic Toolbox

MATLAB’s symbolic toolbox provides a completely separate computer algebra system called Mupad which can handle nonlinear least squares fitting via its stats::reg function.  Here’s how to solve our problem in this environment.

First, you need to start up the mupad application.  You can do this by entering

into MATLAB. Once you are in mupad, the code looks like this

xdata := [-2,-1.64,-1.33,-0.7,0,0.45,1.2,1.64,2.32,2.9]:
ydata := [0.699369,0.700462,0.695354,1.03905,1.97389,2.41143,1.91091,0.919576,-0.730975,-1.42001]:
stats::reg(xdata,ydata,p1*cos(p2*x)+p2*sin(p1*x),[x],[p1,p2],StartingValues=[1,0.2])

The result returned is

[[1.88185085, 0.7002298172], 0.05381269642]

These are our fitted parameters,p1 and p2, along with the sum of squared residuals. The documentation tells us that the optimisation algorithm is Levenberg-Marquardt– this is rather better than the simplex algorithm used by basic MATLAB’s fminsearch.

MATLAB with the NAG Toolbox

The NAG Toolbox for MATLAB is a commercial product offered by the UK based Numerical Algorithms Group.  Their main products are their C and Fortran libraries but they also have a comprehensive MATLAB toolbox that contains something like 1500+ functions.  My University has a site license for pretty much everything they put out and we make great use of it all.  One of the benefits of the NAG toolbox over those offered by The Mathworks is speed.  NAG is often (but not always) faster since its based on highly optimized, compiled Fortran.  One of the problems with the NAG toolbox is that it is difficult to use compared to Mathworks toolboxes.

In an earlier blog post, I showed how to create wrappers for the NAG toolbox to create an easy to use interface for basic nonlinear curve fitting.  Here’s how to solve our problem using those wrappers.

format long
format compact
xdata = [-2,-1.64,-1.33,-0.7,0,0.45,1.2,1.64,2.32,2.9];
ydata = [0.699369,0.700462,0.695354,1.03905,1.97389,2.41143,1.91091,0.919576,-0.730975,-1.42001];

%Note that we don't have to explicitly calculate residuals
fun = @(p,xdata) p(1)*cos(p(2)*xdata)+p(2)*sin(p(1)*xdata);
start = [1,0.2];
[p,fminres]=nag_lsqcurvefit(fun,start,xdata,ydata)

which gives

Warning: nag_opt_lsq_uncon_mod_func_easy (e04fy) returned a warning indicator (5)
p =
1.881850904268710   0.700229557886739
fminres =
0.053812696425390

For convenience, here’s the two files you’ll need to run the above (you’ll also need the NAG Toolbox for MATLAB of course)

MATLAB with curve fitting toolbox

One would expect the curve fitting toolbox to be able to fit such a simple curve and one would be right :)

xdata = [-2,-1.64,-1.33,-0.7,0,0.45,1.2,1.64,2.32,2.9];
ydata = [0.699369,0.700462,0.695354,1.03905,1.97389,2.41143,1.91091,0.919576,-0.730975,-1.42001];

opt = fitoptions('Method','NonlinearLeastSquares',...
'Startpoint',[1,0.2]);
f = fittype('p1*cos(p2*x)+p2*sin(p1*x)','options',opt);

fitobject = fit(xdata',ydata',f)

Note that, unlike every other Mathworks method shown here, xdata and ydata have to be column vectors. The result looks like this

fitobject =

General model:
fitobject(x) = p1*cos(p2*x)+p2*sin(p1*x)
Coefficients (with 95% confidence bounds):
p1 =       1.882  (1.819, 1.945)
p2 =      0.7002  (0.6791, 0.7213)

fitobject is of type cfit:

class(fitobject)

ans =
cfit

In this case it contains two fields, p1 and p2, which are the parameters we are looking for

>> fieldnames(fitobject)
ans =
'p1'
'p2'
>> fitobject.p1
ans =
1.881848414551983
>> fitobject.p2
ans =
0.700229137656802

For maximum information, call the fit command like this:

[fitobject,gof,output] = fit(xdata',ydata',f)

fitobject =

General model:
fitobject(x) = p1*cos(p2*x)+p2*sin(p1*x)
Coefficients (with 95% confidence bounds):
p1 =       1.882  (1.819, 1.945)
p2 =      0.7002  (0.6791, 0.7213)

gof =

sse: 0.053812696487326
rsquare: 0.995722238905101
dfe: 8
rmse: 0.082015773244637

output =

numobs: 10
numparam: 2
residuals: [10x1 double]
Jacobian: [10x2 double]
exitflag: 3
firstorderopt: 3.582047395989108e-05
iterations: 6
funcCount: 21
cgiterations: 0
algorithm: 'trust-region-reflective'
message: [1x86 char]
1. For completeness, you missed the |lsqnonlin| function from the Statistics toolbox. It is basically equivalent to |lsqcurvefit| for solving nonlinear LSQ problems, only exposing a slightly different interface.

% objective function
fun = @(p) (p(1)*cos(p(2)*xdata)+p(2)*sin(p(1)*xdata)) – ydata;

% optimization options
opts = optimset(‘Algorithm’,’levenberg-marquardt’, ‘Display’,’iter’);

% fit
[p,~,res,~,info] = lsqnonlin(fun, pguess, [], [], opts)
fminres = sum(res.^2)

2. Also the Statistics Toolbox has a GUI tool to interactively show the fit and confidence intervals (you can get all that info programmatically as well with nlinfit/nlparci/nlpredci).

screenshot: http://i.imgur.com/ZkTvLfn.png

modelFcn = @(p,xdata) (p(1)*cos(p(2)*xdata)+p(2)*sin(p(1)*xdata));
nlintool(xdata, ydata, modelFcn, pguess)

3. @Amro – Thanks for the mention of lsqnonlin…a usefeul extra interface. I think, however, that lsqnonlin is in the optimisation toolbox – http://www.mathworks.co.uk/help/optim/ug/lsqnonlin.html

4. Hi
The topic and your explanations are fantastic. Would you please guide me to have a good comparsion about accurancy and precision of above methods? I have to choose just one of them to solve my problem and I have to choose the best one.

Thanks a lot.
Mani

5. Hi

I’m beginner in MATLAB but above contents are so understandable for me. I have some experimental points and I used nlinfit function to solve my problem and find the best LS curve fitting result. When I use the nlinfit, it pass MSE= 1.5873e+07 and the result is not very good. I pass my data as below, please check output of nlinfit function to see the output.

X Y
6.1131526e+02 5.7158041e+03
9.1697289e+02 3.3833625e+03
1.2226305e+03 5.2071590e+03
1.5282882e+03 1.3882140e+04
1.8339458e+03 1.4588882e+04
2.1396034e+03 1.2574467e+04
2.4452610e+03 1.0234187e+04
2.7509187e+03 1.1699344e+04
3.0565763e+03 1.4908923e+04
3.3622339e+03 1.1326082e+04
3.6678916e+03 1.2137376e+04
3.9735492e+03 1.1410883e+04
4.2792068e+03 8.2865069e+03
4.5848645e+03 8.7873502e+03
4.8905221e+03 1.0569319e+04
5.1961797e+03 1.0313474e+04
5.5018374e+03 1.1494689e+04
5.8074950e+03 1.0843056e+04
6.1131526e+03 1.0750155e+04
6.4188102e+03 1.4364372e+04
6.7244679e+03 1.1821638e+04
7.0301255e+03 1.0629562e+04
7.3357831e+03 1.0471315e+04
7.6414408e+03 1.1587416e+04
7.9470984e+03 1.2606092e+04
8.2527560e+03 1.0150166e+04
8.5584137e+03 1.1029418e+04
8.8640713e+03 1.4220133e+04
9.1697289e+03 1.6532243e+04
9.4753865e+03 1.4887368e+04
9.7810442e+03 1.3625925e+04
1.0086702e+04 1.4087955e+04
1.0392359e+04 1.3402600e+04
1.0698017e+04 9.5341747e+03
1.1003675e+04 7.8466496e+03
1.1309332e+04 1.2912466e+04
1.1614990e+04 1.6941667e+04
1.1920648e+04 1.4799435e+04
1.2226305e+04 1.2276159e+04
1.2531963e+04 1.1704036e+04
1.2837620e+04 1.1169688e+04
1.3143278e+04 9.5802194e+03
1.3448936e+04 1.3855701e+04
1.3754593e+04 1.4625337e+04
1.4060251e+04 1.1649321e+04
1.4365909e+04 1.3083211e+04
1.4671566e+04 1.4033966e+04
1.4977224e+04 1.6632241e+04
1.5282882e+04 1.4784819e+04
1.5588539e+04 1.4262199e+04
1.5894197e+04 1.5194425e+04
1.6199854e+04 1.4752997e+04
1.6505512e+04 1.2789902e+04
1.6811170e+04 1.3508096e+04
1.7116827e+04 1.6427782e+04
1.7422485e+04 1.6708900e+04
1.7728143e+04 1.7035697e+04
1.8033800e+04 1.5987009e+04
1.8339458e+04 1.7314476e+04
1.8645115e+04 1.7209297e+04
1.8950773e+04 1.5575316e+04
1.9256431e+04 1.5868321e+04
1.9562088e+04 1.5010388e+04
1.9867746e+04 1.4256178e+04
2.0173404e+04 1.3205388e+04
2.0479061e+04 1.3827122e+04
2.0784719e+04 1.5906597e+04
2.1090377e+04 1.8500192e+04
2.1396034e+04 1.8377195e+04
2.1701692e+04 1.6317988e+04
2.2007349e+04 1.5948360e+04
2.2313007e+04 1.9093461e+04
2.2618665e+04 1.8634504e+04
2.2924322e+04 1.4497356e+04
2.3229980e+04 1.5152206e+04
2.3535638e+04 1.8664801e+04
2.3841295e+04 1.9539744e+04
2.4146953e+04 1.6633511e+04
2.4452610e+04 1.8243771e+04
2.4758268e+04 1.8584322e+04
2.5063926e+04 1.6043528e+04
2.5369583e+04 1.8269886e+04
2.5675241e+04 1.9681765e+04
2.5980899e+04 1.5303843e+04
2.6286556e+04 1.6466968e+04
2.6592214e+04 2.0627977e+04
2.6897871e+04 2.0493475e+04
2.7203529e+04 1.5797242e+04
2.7509187e+04 1.3051725e+04

How I should set function’s arguments to reach the best results? Please help me to solve my problem.

Best regards
Parsebay

6. Hi Mike,
I’m experiencing some troubles with curve fitting and I just came across your site, which really makes it way clearer to me but the data I am working with is a bit troublesome and I am not able to find a way to make it work…

let’s say
my Y data is a set of numeric measurements I made: 1×10 vector
my X data is a set of images 1024×1024 px: vector of 1×10 Images
my function F(X,p1,p2) is the central value of the matrix that comes out of the DEconvolution of each image x(1)..x(10), with a Kernel K1,
where K1 is also a 1024×1024 px matrix such that K1=p(1)*exp(-p(2)*Rij)/Rij^2, and Rij is the distance of a pixel ij from the centre of the image. (K1 needs to be =0 for Rij=0 due to indetermination).

I have tried by defining K1 separately and then setting my function F as: mean2(ifft2(fft2(X)./fft2(K1))) but then I do not include the parameters p(1), p(2) anywhere in my function and it does not work.

Do you have any idea about how can I make it?