lsqmultinonlin

Global parametric nonlinear regression with shared parameters

lsqmultinonlin Nonlinear least-squares global regression of multiple data sets

A wrapper function for lsqnonlin which allows simultaneous global curve fitting operations performed on multiple datasets using a parametrized model and tries to estimate the parameters of the model. . Because datasets remain distinct, they may or may not "share" parameter values during the fit process. When a parameter is shared, a single parameter value is calculated for all datasets.

lsqmultinonlin attempts to solve problems of the form:
min sum {FUN(X,beta)-Y}.^2 in a least square sense with respect to beta where X and the values returned by FUN can be vectors.

In this model FUN=[f1(X,beta);f2(X,beta);...]

Note that the convergence might be slow due to the simultaneous minimizations

Beta = lsqmultinonlin(FUN,Beta0,X,Y) starts at the vector Beta0 and finds a minimum Beta that fits the functions in FUN to the datasets Y.
FUN accepts the inputs X and returns a vector of function values F evaluated at X. NOTE: FUN should return FUN(X) and not the sum-of-squares sum(FUN(X).^2)). (FUN(X) is summed and squared implicitly in the algorithm.)

INPUT:
X,Y: Matrices containing the X,Y vectors of the fitted data sets.
X=[X1, X2, ...], Y=[Y1, Y2, ...]

FUN: Cell array containing model functions for each data set as functions of the column vectors X1, X2, ...

beta0: Vector containing initial guess of the fitting parameters.

OUTPUT:
Beta: set of shared or unshared parameters the minimize the error

Beta = lsqmultinonlin(FUN,Beta0,X,Y,LB,UB) defines a set of lower and upper bounds on the design parameters, Beta, so that the solution is in the range LB <= Beta <= UB. Use empty matrices for LB and UB if no bounds exist. Set LB(i)= -Inf if Beta(i) is unbounded below; set
UB(i) = Inf if Beta(i) isunbounded above.

Beta = lsqmultinonlin(FUN,X0,X,Y,LB,UB,OPTIONS) minimizes with the default optimization parameters replaced by values in OPTIONS, an argument created with the OPTIMOPTIONS function. See OPTIMOPTIONS for details.

[Beta,RESNORM] = lsqmultinonlin(FUN,Beta0,X,Y,...) returns the value of the squared 2-norm of the residual at Beta: sum(FUN(Beta,X)-Y).^2.

[Beta,RESNORM,RESIDUAL] = lsqmultinonlin(FUN,Beta0,X,Y,...) returns the value of the residual at the solution Beta: RESIDUAL = FUN(Beta,X)-Y.

[Beta,RESNORM,RESIDUAL,EXITFLAG] = lsqmultinonlin(FUN,Beta0,X,Y,...) returns an EXITFLAG that describes the exit condition. Possible values of EXITFLAG and the corresponding exit conditions are listed below. See the documentation of lsqnonlin for a complete description.

1 lsqnonlin converged to a solution.
2 Change in X too small.
3 Change in RESNORM too small.
4 Computed search direction too small.
0 Too many function evaluations or iterations.
-1 Stopped by output/plot function.
-2 Bounds are inconsistent.

[Beta,RESNORM,RESIDUAL,EXITFLAG,OUTPUT] = lsqmultinonlin(FUN,Beta0,X,Y,...) returns a structure OUTPUT with the number of iterations taken in OUTPUT.iterations, the number of function evaluations in OUTPUT.funcCount, the algorithm used in OUTPUT.algorithm, the number of CG iterations (if used) in OUTPUT.cgiterations, the first-order optimality (if used) in OUTPUT.firstorderopt, and the exit message in OUTPUT.message.

[Beta,RESNORM,RESIDUAL,EXITFLAG,OUTPUT,LAMBDA] = lsqmultinonlin(FUN,Beta0,X,Y,...) returns the set of Lagrangian multipliers, LAMBDA, at the solution: LAMBDA.lower for LB and LAMBDA.upper for UB.

[Beta,RESNORM,RESIDUAL,EXITFLAG,OUTPUT,LAMBDA,JACOBIAN] = lsqmultinonlin(FUN,Beta0,X,Y,...) returns the Jacobian of FUN-Y at Beta.

EXAMPLE:
% Generate X vectors for both data sets
X1 = 0:1:100;
X2 = 0:1:100;
X=[X1',X2'];
% Generate Y data with some noise
Y1 = cos(2*pi*0.5*X1).*exp(-X2/5) + 0.05*randn(size(X1));
Y2 = 0.5 + 2*exp(-(X2/5)) + 0.05*randn(size(X2));
Y=[Y1',Y2'];
% Define fitting functions and parameters, with identical
% exponential decay for both data sets
mdl1 = @(Beta,X) cos(2*pi*Beta(1)*X(:,1)).*exp(-X(:,2)/Beta(2));
mdl2 = @(Beta,X) Beta(4) + Beta(3)*exp(-(X(:,2)/Beta(2)));
%
% Prepare input for lsqmultinonlin and perform fitting
FUN = {mdl1, mdl2};
Beta0 = [0.2, 1, 1, 1];
[Beta,RESNORM,RESIDUAL,EXITFLAG,OUTPUT,LAMBDA,JACOBIAN] = lsqmultinonlin(FUN, Beta0, X, Y,[0,0,0,0],[1,10,4,1]);
% Plot results
figure;
hold all;
box on;
scatter(X1,Y1);
scatter(X2,Y2);
plot(X1,mdl1(Beta,X),'Color','blue');
plot(X2,mdl2(Beta,X),'Color',[0 0.5 0]);