Nonlinear fit with constraints in R2012b

Hi

I am running R2012b and doing nonlinear weighted regression as explained here: http://www.mathworks.se/help/stats/examples/weighted-nonlinear-regression.html. In other words, I am using NonLinearModel.fit.

However is there a way to constrain one of the fit parameters to be positive?

Best, Niles.

Answer 1

[Star Strider]

The Statistics Toolbox functions (such as you referenced in the link) won't let you constrain the parameters but the Optimization Toolbox functions (such as lsqcurvefit) will.

You can do weighted nonlinear regression with lsqcurvefit, using an interim anonymous function. So if the model you have coded as a function and want to fit is called ‘nonlinmodel’, you would normally call it from lsqcurvefit this way:

x = lsqcurvefit(@nonlinmodel,x0,xdata,ydata,lb,ub);

To do weighted nonlinear regression with it, you need to incorporate your weight vector into an anonymous function, and then call the anonymous fiunction from lsqcurvefit:

Weights = 1./(N.*SE.^2);

nonlinmodelW = @(B,t) Weights .* nonlinearmodel(B,t);

x = lsqcurvefit(nonlinmodelW,x0,xdata,ydata,lb,ub);

For ‘Weights’, I used the typical standard errors usually provided, with N being the number of observations used to calculate each of them. This produces inverse-variance weights, the conventional scheme. Note that ‘Weights’ is a vector of the same dimensions as the vector ‘nonlinearmodel’ returns.

When you simulate your function with your estimated parameters to plot with your data, remember to use the anonymous function. It (here ‘nonlinmodelW’) is the function you fitted, not your original ‘nonlinmodel’.

This anonymous function approach will work with both lsqcurvefit and nlinfit. Before weights were added to nlinfit in 2012b, this was the standard way of doing weighted nonlinear regression with both. I've done weighted nonlinear regression using this approach many times with lsqcurvefit in order to take advantage of the parameter constraints.

[Niles Martinsen]

Hi Star Strider

Thanks for the help. I need uncertainties on the fit parameters too, but as far as I could tell (<http://www.mathworks.se/support/solutions/en/data/1-18QY1/index.html?product=ST&solution=1-18QY1)>, LSCURVEFIT does not support this. Can I work around this somehow?

[Star Strider]

I couldn't find my references for it, so I went online looking and found this enlightening thread from the Newsgroup. He was using lsqnonlin, but lsqcurvefit is the same routine, set up for curve-fitting. I haven't coded his routine yet, but I will.

I extracted his pseudocode snippet here:

%===============================
v=residuals from nonlinear least squares;
J=jacobian;

[rows,cols]=size(J);

% A posteriori variance factor
 Sigma_o = sqrt((v'*v) / (rows-cols));

% % Unscaled Covariance Matrix (Inverse of the Normal Equation matrix)
 Q_xx = inv(J'*J);

% Precision measures ,i.e. the standard deviation of each parameter
 diagonals_of_Q_xx = diag(Q_xx);
Precision_of_solved_parameters = Sigma_o .* sqrt(diagonals_of_Q_xx);
%================================

It seems straightforward. I'll code it and then compare it with the same estimates for nlparci.

Since you don't have the Statistics Toolbox, you might also find these helpful for computing the normal probability from a z-statistic and the z-statistic (or critical value) from a probability:

   P = @(z) erfc(-z/sqrt(2))/2;
   CV = @(alpha) -sqrt(2) * erfcinv(2*alpha/2);

‘CV’ (‘Critical Value’) would be necessary to determine the 95% (or other) confidence intervals from the parameter standard deviations, setting alpha at 0.975 or as a vector [0.025 0.975].

[Zhaorong Wang]

Thanks Star Strider for providing a nice trick to do weighted fitting using the most powerful lsqcurvefit! However, one thing has to be corrected, if you weight the fitting function "nonlinmodel", you have to weight the ydata as well. Because you are minimizing

sum((nonlinmodelW(B,xdata) - ydata).^2),

readers better deal with this carefully, like writing some wrapper function lsqcurvefit_weighted() for lsqcurvefit()

[Sarah Garre]

Dear Zhaorong Wang, thanks for putting our attention to this important detail. How do you suggest to make a 'wrapper'?

Answer 2

[Sarah Garre]

if I understand it well, it should be enough to manage it like this:

Weights = 1./(N.*SE.^2);

nonlinmodelW = @(B,t) Weights .* nonlinearmodel(B,t);

x = lsqcurvefit(nonlinmodelW,x0,xdata,Weights .* ydata,lb,ub);

Am I right?

[Richard Andersson]

Thanks everyone for your insight on this problem. I would just ask that you clarify exactly what the inputs "B" and "t" are.

[Star Strider]

Sarah — That looks correct. I haven’t tested it to be sure.

Richeard — The argument B is the vector of parameters (usually termed β in the literature), and t is the independent variable vector, probably time in this example.

The same objective function format works for both lsqcurvefit and nlinfit, but the arguments must be in the (B,t) order.

In the last few versions, nlinfit has the option of a weighted regression. I haven’t needed to use it recently, so I refer you to the documentation.

Answer 3

[Bhavin Khatri]

I think by usual maximum likelihood arguments the weights should be

weights = 1./se.^2

not

weights 1./(N.*se.^2).

Intuitively, it doesn't make sense either, since according to the latter, the more observations you have for a particular data point the smaller the weight!

(Of course you can choose whatever weights you want, but assuming Gaussian errors etc..., if you want the maximum likelihood estimate the weights should be the inverse of the standard error.)