This example shows how to fit a function to data using
lsqcurvefit together with
Many fitting problems have multiple local solutions.
MultiStart can help find the global solution, meaning the best fit. While you can use
lsqnonlin as the local solver, this example uses
lsqcurvefit simply because it has a convenient syntax.
The model is
where the input data is , and the parameters , , , and are the unknown model coefficients.
Write an anonymous function that takes a data matrix
N rows and two columns, and returns a response vector with
N rows. It also takes a coefficient matrix
p, corresponding to the coefficient vector
fitfcn = @(p,xdata)p(1) + p(2)*xdata(:,1).*sin(p(3)*xdata(:,2)+p(4));
Create 200 data points and responses. Use the values . Include random noise in the response.
rng default % for reproducibility N = 200; % number of data points preal = [-3,1/4,1/2,1]; % real coefficients xdata = 5*rand(N,2); % data points ydata = fitfcn(preal,xdata) + 0.1*randn(N,1); % response data with noise
Set bounds for
lsqcurvefit. There is no reason for
in absolute value, because the sine function takes values in its full range over any interval of width
. Assume that the coefficient
must be smaller than 20 in absolute value, because allowing a very high frequency can cause unstable responses or spurious convergence.
lb = [-Inf,-Inf,-20,-pi]; ub = [Inf,Inf,20,pi];
Set the initial point arbitrarily to (5,5,5,0).
p0 = 5*ones(1,4); % Arbitrary initial point p0(4) = 0; % so the initial point satisfies the bounds
Fit the parameters to the data, starting at
[xfitted,errorfitted] = lsqcurvefit(fitfcn,p0,xdata,ydata,lb,ub)
Local minimum possible. lsqcurvefit stopped because the final change in the sum of squares relative to its initial value is less than the default value of the function tolerance. xfitted = -2.6149 -0.0238 6.0191 -1.6998 errorfitted = 28.2524
lsqcurvefit found a local solution that is not particularly close to the model parameter values (-3,1/4,1/2,1).
Create a problem structure so
MultiStart can solve the same problem.
problem = createOptimProblem('lsqcurvefit','x0',p0,'objective',fitfcn,... 'lb',lb,'ub',ub,'xdata',xdata,'ydata',ydata);
Solve the fitting problem using
MultiStart with 50 iterations. Plot the smallest error as the number of
ms = MultiStart('PlotFcns',@gsplotbestf); [xmulti,errormulti] = run(ms,problem,50)
MultiStart completed the runs from all start points. All 50 local solver runs converged with a positive local solver exit flag. xmulti = -2.9852 -0.2472 -0.4968 -1.0438 errormulti = 1.6464
MultiStart found a global solution near the parameter values
. (This is equivalent to a solution near
, because changing the sign of all the coefficients except the first gives the same numerical values of
fitfcn.) The norm of the residual error decreased from about 28 to about 1.6, a decrease of more than a factor of 10.