Can I use the Jacobian provided by 'lsqnonlin' to compute the confidence intervals using 'nlparci'?
34 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
Cesar de Araujo Filho
il 5 Mag 2019
Commentato: Adam Danz
il 6 Mag 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 Commenti
Risposta accettata
Adam Danz
il 5 Mag 2019
Modificato: Adam Danz
il 6 Mag 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 Commenti
Adam Danz
il 6 Mag 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).
Vedere anche
Categorie
Scopri di più su Interpolation in Help Center e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!