Can I use the Jacobian provided by 'lsqnonlin' to compute the confidence intervals using 'nlparci'?
34 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Cesar de Araujo Filho
el 5 de Mayo de 2019
Comentada: Adam Danz
el 6 de Mayo de 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 comentarios
Respuesta aceptada
Adam Danz
el 5 de Mayo de 2019
Editada: Adam Danz
el 6 de Mayo de 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 comentarios
Adam Danz
el 6 de Mayo de 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).
Ver también
Categorías
Más información sobre Interpolation en Help Center y File Exchange.
Productos
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!