{"id":5196,"date":"2013-12-06T13:23:53","date_gmt":"2013-12-06T12:23:53","guid":{"rendered":"http:\/\/www.walkingrandomly.com\/?p=5196"},"modified":"2013-12-06T16:13:35","modified_gmt":"2013-12-06T15:13:35","slug":"simple-nonlinear-least-squares-curve-fitting-in-matlab","status":"publish","type":"post","link":"https:\/\/walkingrandomly.com\/?p=5196","title":{"rendered":"Simple nonlinear least squares curve fitting in MATLAB"},"content":{"rendered":"<p>A question I get asked a lot is &#8216;How can I do <a href=\"http:\/\/en.wikipedia.org\/wiki\/Non-linear_least_squares\">nonlinear least squares<\/a> curve fitting in X?&#8217; where X might be MATLAB, Mathematica or a whole host of alternatives. \u00a0Since this is such a common query, I thought I&#8217;d write up how to do it for a very simple problem in several systems that I&#8217;m interested in<\/p>\n<p>This is the MATLAB version. For other versions,see the list below<\/p>\n<ul>\n<li><a href=\"https:\/\/www.walkingrandomly.com\/?p=5181\">Simple nonlinear least squares curve fitting in Julia<\/a><\/li>\n<li><a href=\"https:\/\/www.walkingrandomly.com\/?p=5218\">Simple nonlinear least squares curve fitting in Maple<\/a><\/li>\n<li><a href=\"https:\/\/www.walkingrandomly.com\/?p=5180\">Simple nonlinear least squares curve fitting in Mathematica<\/a><\/li>\n<li><a href=\"https:\/\/www.walkingrandomly.com\/?p=5215\">Simple nonlinear least squares curve fitting in Python<\/a><\/li>\n<li><a href=\"https:\/\/www.walkingrandomly.com\/?p=5254\">Simple nonlinear least squares curve fitting in R<\/a><\/li>\n<\/ul>\n<p><strong>The problem<\/strong><\/p>\n<pre>xdata = -2,-1.64,-1.33,-0.7,0,0.45,1.2,1.64,2.32,2.9\r\nydata = 0.699369,0.700462,0.695354,1.03905,1.97389,2.41143,1.91091,0.919576,-0.730975,-1.42001<\/pre>\n<p>and you&#8217;d like to fit the function<\/p>\n<p><img loading=\"lazy\" decoding=\"async\" src=\"https:\/\/walkingrandomly.com\/wp-content\/ql-cache\/quicklatex.com-17f46b025dd91e024cf2dc04211e19ba_l3.png\" class=\"ql-img-inline-formula quicklatex-auto-format\" alt=\"&#32;&#70;&#40;&#112;&#95;&#49;&#44;&#112;&#95;&#50;&#44;&#120;&#41;&#32;&#61;&#32;&#112;&#95;&#49;&#32;&#92;&#99;&#111;&#115;&#40;&#112;&#95;&#50;&#32;&#120;&#41;&#43;&#112;&#95;&#50;&#32;&#92;&#115;&#105;&#110;&#40;&#112;&#95;&#49;&#32;&#120;&#41;&#32;\" title=\"Rendered by QuickLaTeX.com\" height=\"18\" width=\"297\" style=\"vertical-align: -4px;\"\/><\/p>\n<p>using nonlinear least squares. \u00a0You&#8217;re starting guesses for the parameters are p1=1 and P2=0.2<\/p>\n<p>For now, we are primarily interested in the following results:<\/p>\n<ul>\n<li>The fit parameters<\/li>\n<li>Sum of squared residuals<\/li>\n<\/ul>\n<p>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.<\/p>\n<p>How you proceed depends on which toolboxes you have. Contrary to popular belief, you don&#8217;t need the Curve Fitting toolbox to do curve fitting&#8230;particularly when the fit in question is as basic as this. Out of the 90+ toolboxes sold by The Mathworks, I&#8217;ve only been able to look through the subset I have access to so I may have missed some alternative solutions.<\/p>\n<p><strong>Pure MATLAB solution (No toolboxes)<\/strong><\/p>\n<p>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. \u00a0Basic MATLAB comes with the fminsearch function which is based on the <a href=\"http:\/\/en.wikipedia.org\/wiki\/Nelder%E2%80%93Mead_method\">Nelder-Mead simplex method<\/a>. \u00a0For this particular problem, it works OK but will not be suitable for more complex fitting problems. \u00a0Here&#8217;s the code<\/p>\n<pre>format compact\r\nformat long\r\nxdata = [-2,-1.64,-1.33,-0.7,0,0.45,1.2,1.64,2.32,2.9];\r\nydata = [0.699369,0.700462,0.695354,1.03905,1.97389,2.41143,1.91091,0.919576,-0.730975,-1.42001];\r\n\r\n%Function to calculate the sum of residuals for a given p1 and p2\r\nfun = @(p) sum((ydata - (p(1)*cos(p(2)*xdata)+p(2)*sin(p(1)*xdata))).^2);\r\n\r\n%starting guess\r\npguess = [1,0.2];\r\n\r\n%optimise\r\n[p,fminres] = fminsearch(fun,pguess)<\/pre>\n<p>This gives the following results<\/p>\n<pre>p =\r\n1.881831115804464 0.700242006994123\r\nfminres =\r\n0.053812720914713<\/pre>\n<p>All we get here are the parameters and the sum of squares of the residuals. \u00a0If you want more information such as 95% confidence intervals, you&#8217;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.<\/p>\n<p><strong>MATLAB with optimisation toolbox<\/strong><\/p>\n<p>With respect to this problem, the optimisation toolbox gives you two main advantages over pure MATLAB. \u00a0The first is that better optimisation routines are available so more complex problems (such as those with constraints) can be solved and in less time. \u00a0The second is the provision of the lsqcurvefit function which is specifically designed to solve curve fitting problems. \u00a0Here&#8217;s the code<\/p>\n<pre>format long\r\nformat compact\r\nxdata = [-2,-1.64,-1.33,-0.7,0,0.45,1.2,1.64,2.32,2.9];\r\nydata = [0.699369,0.700462,0.695354,1.03905,1.97389,2.41143,1.91091,0.919576,-0.730975,-1.42001];\r\n\r\n%Note that we don't have to explicitly calculate residuals\r\nfun = @(p,xdata) p(1)*cos(p(2)*xdata)+p(2)*sin(p(1)*xdata);\r\n\r\n%starting guess\r\npguess = [1,0.2];\r\n\r\n[p,fminres] = lsqcurvefit(fun,pguess,xdata,ydata)<\/pre>\n<p>This gives the results<\/p>\n<pre>p =\r\n   1.881848414551983   0.700229137656802\r\nfminres =\r\n   0.053812696487326<\/pre>\n<p><strong>MATLAB with statistics toolbox<\/strong><\/p>\n<p>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<\/p>\n<pre>%set up for both fit commands in the stats toolbox\r\nxdata = [-2,-1.64,-1.33,-0.7,0,0.45,1.2,1.64,2.32,2.9];\r\nydata = [0.699369,0.700462,0.695354,1.03905,1.97389,2.41143,1.91091,0.919576,-0.730975,-1.42001];\r\n\r\nfun = @(p,xdata) p(1)*cos(p(2)*xdata)+p(2)*sin(p(1)*xdata);\r\n\r\npguess = [1,0.2];<\/pre>\n<p>Method 1 makes use of <strong>NonLinearModel.fit<\/strong><\/p>\n<pre>mdl = NonLinearModel.fit(xdata,ydata,fun,pguess)<\/pre>\n<p>The returned object is a <strong>NonLinearModel<\/strong> object<\/p>\n<pre>class(mdl)\r\nans =\r\nNonLinearModel<\/pre>\n<p>which contains all sorts of useful stuff<\/p>\n<pre>mdl = \r\nNonlinear regression model:\r\n    y ~ p1*cos(p2*xdata) + p2*sin(p1*xdata)\r\n\r\nEstimated Coefficients:\r\n          Estimate             SE                 \r\n    p1      1.8818508110535      0.027430139389359\r\n    p2    0.700229815076442    0.00915260662357553\r\n\r\n          tStat               pValue              \r\n    p1    68.6052223191956    2.26832562501304e-12\r\n    p2    76.5060538352836    9.49546284187105e-13\r\n\r\nNumber of observations: 10, Error degrees of freedom: 8\r\nRoot Mean Squared Error: 0.082\r\nR-Squared: 0.996,  Adjusted R-Squared 0.995\r\nF-statistic vs. zero model: 1.43e+03, p-value = 6.04e-11<\/pre>\n<p>If you don&#8217;t need such heavyweight infrastructure, you can make use of the statistic toolbox&#8217;s <strong>nlinfit<\/strong> function<\/p>\n<pre>[p,R,J,CovB,MSE,ErrorModelInfo] = nlinfit(xdata,ydata,fun,pguess);<\/pre>\n<p>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).<\/p>\n<p>Both <strong>nlinfit<\/strong> and <strong>NonLinearModel.fit<\/strong> use the same minimisation algorithm which is based on Levenberg-Marquardt<\/p>\n<p><strong>MATLAB with Symbolic Toolbox<\/strong><\/p>\n<p>MATLAB&#8217;s symbolic toolbox provides a completely separate computer algebra system called Mupad which can handle nonlinear least squares fitting via its stats::reg function. \u00a0Here&#8217;s how to solve our problem in this environment.<\/p>\n<p>First, you need to start up the mupad application. \u00a0You can do this by entering<\/p>\n<pre>mupad<\/pre>\n<p>into MATLAB. Once you are in mupad, the code looks like this<\/p>\n<pre>xdata := [-2,-1.64,-1.33,-0.7,0,0.45,1.2,1.64,2.32,2.9]:\r\nydata := [0.699369,0.700462,0.695354,1.03905,1.97389,2.41143,1.91091,0.919576,-0.730975,-1.42001]:\r\nstats::reg(xdata,ydata,p1*cos(p2*x)+p2*sin(p1*x),[x],[p1,p2],StartingValues=[1,0.2])<\/pre>\n<p>The result returned is<\/p>\n<pre>[[1.88185085, 0.7002298172], 0.05381269642]<\/pre>\n<p>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&#8211; this is rather better than the simplex algorithm used by basic MATLAB&#8217;s fminsearch.<\/p>\n<p><strong>MATLAB with the NAG Toolbox<\/strong><\/p>\n<p>The <a href=\"http:\/\/www.nag.co.uk\/numeric\/MB\/start.asp\">NAG Toolbox for MATLAB<\/a> is a commercial product offered by the UK based Numerical Algorithms Group. \u00a0Their main products are their <a href=\"http:\/\/www.nag.co.uk\/numeric\/CL\/CLdescription.asp\">C<\/a> and <a href=\"http:\/\/www.nag.co.uk\/numeric\/fl\/FLdescription.asp\">Fortran<\/a> libraries but they also have a comprehensive MATLAB toolbox that contains something like 1500+ functions. \u00a0My University has a site license for pretty much everything they put out and we make great use of it all. \u00a0One of the benefits of the NAG toolbox over those offered by The Mathworks is speed. \u00a0NAG is often (but not always) faster since its based on highly optimized, compiled Fortran. \u00a0One of the problems with the NAG toolbox is that it is difficult to use compared to Mathworks toolboxes.<\/p>\n<p>In an earlier blog post, I <a href=\"https:\/\/www.walkingrandomly.com\/?p=3535\">showed how to create wrappers for the NAG toolbox<\/a> to create an easy to use interface for basic nonlinear curve fitting. \u00a0Here&#8217;s how to solve our problem using those wrappers.<\/p>\n<pre>format long\r\nformat compact\r\nxdata = [-2,-1.64,-1.33,-0.7,0,0.45,1.2,1.64,2.32,2.9];\r\nydata = [0.699369,0.700462,0.695354,1.03905,1.97389,2.41143,1.91091,0.919576,-0.730975,-1.42001];\r\n\r\n%Note that we don't have to explicitly calculate residuals\r\nfun = @(p,xdata) p(1)*cos(p(2)*xdata)+p(2)*sin(p(1)*xdata);\r\nstart = [1,0.2];\r\n[p,fminres]=nag_lsqcurvefit(fun,start,xdata,ydata)<\/pre>\n<p>which gives<\/p>\n<pre>Warning: nag_opt_lsq_uncon_mod_func_easy (e04fy) returned a warning indicator (5) \r\np =\r\n   1.881850904268710   0.700229557886739\r\nfminres =\r\n   0.053812696425390<\/pre>\n<p>For convenience, here&#8217;s the two files you&#8217;ll need to run the above (you&#8217;ll also need the NAG Toolbox for MATLAB of course)<\/p>\n<ul>\n<li><a href=\"https:\/\/www.walkingrandomly.com\/images\/NAG\/MATLAB\/lsqcurvefit\/1\/nag_lsqcurvefit.m\">nag_lsqcurvefit.m<\/a><\/li>\n<li><a href=\"https:\/\/www.walkingrandomly.com\/images\/NAG\/MATLAB\/lsqcurvefit\/1\/nag_lsqcurvefit_aux.m\">nag_lsqcurvefit_aux.m<\/a><\/li>\n<\/ul>\n<p><strong>MATLAB with curve fitting toolbox<\/strong><\/p>\n<p>One would expect the curve fitting toolbox to be able to fit such a simple curve and one would be right :)<\/p>\n<pre>xdata = [-2,-1.64,-1.33,-0.7,0,0.45,1.2,1.64,2.32,2.9];\r\nydata = [0.699369,0.700462,0.695354,1.03905,1.97389,2.41143,1.91091,0.919576,-0.730975,-1.42001];\r\n\r\nopt = fitoptions('Method','NonlinearLeastSquares',...\r\n               'Startpoint',[1,0.2]);\r\nf = fittype('p1*cos(p2*x)+p2*sin(p1*x)','options',opt);\r\n\r\nfitobject = fit(xdata',ydata',f)<\/pre>\n<p>Note that, unlike every other Mathworks method shown here, xdata and ydata have to be column vectors. The result looks like this<\/p>\n<pre>fitobject = \r\n\r\n     General model:\r\n     fitobject(x) = p1*cos(p2*x)+p2*sin(p1*x)\r\n     Coefficients (with 95% confidence bounds):\r\n       p1 =       1.882  (1.819, 1.945)\r\n       p2 =      0.7002  (0.6791, 0.7213)<\/pre>\n<p>fitobject is of type <a href=\"http:\/\/www.mathworks.co.uk\/help\/curvefit\/cfit.html\">cfit<\/a>:<\/p>\n<pre> class(fitobject)\r\n\r\nans =\r\ncfit<\/pre>\n<p>In this case it contains two fields, p1 and p2, which are the parameters we are looking for<\/p>\n<pre>&gt;&gt; fieldnames(fitobject)\r\nans = \r\n    'p1'\r\n    'p2'\r\n&gt;&gt; fitobject.p1\r\nans =\r\n   1.881848414551983\r\n&gt;&gt; fitobject.p2\r\nans =\r\n    0.700229137656802<\/pre>\n<p>For maximum information, call the fit command like this:<\/p>\n<pre>[fitobject,gof,output] = fit(xdata',ydata',f)\r\n\r\nfitobject = \r\n\r\n     General model:\r\n     fitobject(x) = p1*cos(p2*x)+p2*sin(p1*x)\r\n     Coefficients (with 95% confidence bounds):\r\n       p1 =       1.882  (1.819, 1.945)\r\n       p2 =      0.7002  (0.6791, 0.7213)\r\n\r\ngof = \r\n\r\n           sse: 0.053812696487326\r\n       rsquare: 0.995722238905101\r\n           dfe: 8\r\n    adjrsquare: 0.995187518768239\r\n          rmse: 0.082015773244637\r\n\r\noutput = \r\n\r\n           numobs: 10\r\n         numparam: 2\r\n        residuals: [10x1 double]\r\n         Jacobian: [10x2 double]\r\n         exitflag: 3\r\n    firstorderopt: 3.582047395989108e-05\r\n       iterations: 6\r\n        funcCount: 21\r\n     cgiterations: 0\r\n        algorithm: 'trust-region-reflective'\r\n          message: [1x86 char]<\/pre>\n","protected":false},"excerpt":{"rendered":"<p>A question I get asked a lot is &#8216;How can I do nonlinear least squares curve fitting in X?&#8217; where X might be MATLAB, Mathematica or a whole host of alternatives. \u00a0Since this is such a common query, I thought I&#8217;d write up how to do it for a very simple problem in several systems [&hellip;]<\/p>\n","protected":false},"author":3,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"jetpack_post_was_ever_published":false,"footnotes":"","jetpack_publicize_message":"","jetpack_publicize_feature_enabled":true,"jetpack_social_post_already_shared":true,"jetpack_social_options":{"image_generator_settings":{"template":"highway","default_image_id":0,"font":"","enabled":false},"version":2}},"categories":[4,11],"tags":[],"class_list":["post-5196","post","type-post","status-publish","format-standard","hentry","category-math-software","category-matlab"],"jetpack_publicize_connections":[],"jetpack_featured_media_url":"","jetpack_shortlink":"https:\/\/wp.me\/p3swhs-1lO","jetpack_sharing_enabled":true,"_links":{"self":[{"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/posts\/5196","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/users\/3"}],"replies":[{"embeddable":true,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=%2Fwp%2Fv2%2Fcomments&post=5196"}],"version-history":[{"count":16,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/posts\/5196\/revisions"}],"predecessor-version":[{"id":5287,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/posts\/5196\/revisions\/5287"}],"wp:attachment":[{"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=%2Fwp%2Fv2%2Fmedia&parent=5196"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=%2Fwp%2Fv2%2Fcategories&post=5196"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=%2Fwp%2Fv2%2Ftags&post=5196"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}