Can I use the Jacobian provided by 'lsqnonlin' to compute the confidence intervals using 'nlparci'?
34 views (last 30 days)
Show older comments
Cesar de Araujo Filho
on 5 May 2019
Commented: Adam Danz
on 6 May 2019
The function 'nlparci' accepts as input the Jacobian of the regressing function at the optimal point. In the explanation of this function, it is only mentioned the Jacobian which is obtained from 'nlinfit'. Nevertheless, 'lsqnonlin' also provides a Jacobian output. Is it correct to use the Jacobian from 'lsqnonlin' directly in 'nlparci'?
0 Comments
Accepted Answer
Adam Danz
on 5 May 2019
Edited: Adam Danz
on 6 May 2019
[estParams,~,residual,~,~,~,jacobian] = lsqnonlin();
CI = nlparci(estParams,residual,'jacobian',jacobian);
An alternative is to calculate the covariance matrix to calculate the confidence intervals directly (shown below). There are typically very small differences in the results between the two methods.
lastwarn(''); %Clear warning memory
alpha = 0.05; %significance level
df = length(residual) - numel(estParams); %degrees of freedom
crit = tinv(1-alpha/2,df); %critical value
covm = inv(jacobian'*jacobian) * var(residual); %covariance matrix
[~, warnId] = lastwarn; %detect if inv threw 'nearly singular' warning.
covmIdx = sub2ind(size(covm),1:size(covm,1),1:size(covm,2)); %indices of the diag of covm
CI = nan(numel(estParams),2);
if ~strcmp(warnId, 'MATLAB:nearlySingularMatrix')
CI(:,1) = estParams - crit * sqrt(covm(covmIdx));
CI(:,2) = estParams + crit * sqrt(covm(covmIdx));
end
2 Comments
Adam Danz
on 6 May 2019
I wrote that code a while back for my own needs and remember comparing the results to nlparci() and finding very miniscule differences but not precisely the same values (differences several orders magnitude smaller than the data). .
The mistake you mentioned reminds me of a recent push I'm making to use numel() instead of length(). You can see I used both in my code (which was written a while back).
See Also
Categories
Find more on Interpolation in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!