Uncertainy in constrained lsqlin fitting

Hi,

I need to perform a linear regression with lower bounds on the parameters (they have to be positive or zero). I used

lb = zeros(1,13); B = lsqlin(matrix,data,[ ],[ ],[ ],[ ],lb,[ ]);

where data=matrix*B.

How can I estimate the uncertainty on the B parameters ? Using the function 'regress' I can have the uncertainty, but no constraints. Using lsqlin I can set constraints, but how to get the uncertainty ?

Many thanks

Paola

Answer 1

[Matt J]

Once you've solved the problem with lsqlin, you can transform the problem to the equivalent unconstrained one

     min.  norm(matrix*Z.^2-data)

where Z.^2=B is a change of variables and then re-solve. You can re-solve with fminunc or another unconstrained solver in zero iterations by initializing with Z=sqrt(B), where B is the solution already obtained with lsqlin. Because this is an unconstrained problem, though, you can now request uncertainty information on Z from the solver. For example, if you use fminunc, you can output the Hessian and use it to calculate uncertainties on Z in the usual ways,

    [x,fval,exitflag,output,grad,hessian] = fminunc(...)

Finally, once you have the uncertainty of Z, you can translate it back to uncertainties on B using the change of variables formula Z.^2=B. So, I suppose the relationship would be B_errors = sqrt(Z_errors).

[Sterling Baird]

Do you have an example of how to use the Hessian "to calculate uncertainties on Z in the usual ways"? I'll plan to try parsing through some of the stack exchange questions on this topic (e.g. take the inverse of the Hessian to get covariance matrix). Any references for this procedure you could point me to? (i.e. the change of variables you mentioned for transforming the problem)

[Matt J]

Do you have an example of how to use the Hessian "to calculate uncertainties on Z in the usual ways"?

See section 28 and 29 of John D'Errico's Optimization Tips and Tricks,

https://www.mathworks.com/matlabcentral/fileexchange/8553-optimization-tips-and-tricks?s_tid=srchtitle

In this case the Hessian would be the estimate of J.'*J.

Any references for this procedure you could point me to? (i.e. the change of variables you mentioned for transforming the problem)

No, not that I'm aware of.