Can I use the Jacobian provided by 'lsqnonlin' to compute the confidence intervals using 'nlparci'?

41 次查看(过去 30 天)
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'?

采纳的回答

Adam Danz
Adam Danz 2019-5-5
编辑:Adam Danz 2019-5-6
Yes, you can use it but see Matt J's warning about using bounded constraints.
[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 个评论
Cesar de Araujo Filho
Thank you Adam! Your code was providential. The reason why I asked this question was because I wanted to write a code that could calculate the confidence interval without the need of the any Toolbox. But the code I wrote did not give me the same results as the 'nlparci' function. Then I started questioning every possible reason for it.
Your code gives exactly the same result as the 'nlparci' but with the advantage that I can see what's happening behind. I now come to realize that my code was not actually "wrong". I had defined the 'residual' matrix NOT as a collumn matrix, but as a rectangular one. In this case, the degrees of freedom estimated by the 'nlparci' were significantly smaller because they were taking into account only one dimension of the 'residual' matrix. Here is how my code goes:
function result=UDnlparci(par,res,J,alfa)
df=numel(res)-numel(par)
s2=(norm(res)).^2/df;
nu=s2*inv(J'*J);
tDist=@(t) gamma((df+1)./2)./gamma(df./2)*1./sqrt(df.*pi()).*(1+t.^2/df).^(-(df+1)./2);
fmin=@(x) -(1-alfa./2)+integral(tDist,-Inf,x,'RelTol',1E-9,'AbsTol',1e-12);
options=optimset('Display','off','tolx',1E-12,'tolfun',1E-12);
tinv=fzero(fmin,1,options);
result(:,1)=par'-tinv*sqrt(diag(nu));
result(:,2)=par'+tinv*sqrt(diag(nu));
end
Thanks a lot for the prompt answer!
Adam Danz
Adam Danz 2019-5-6
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).

请先登录,再进行评论。

更多回答(1 个)

Matt J
Matt J 2019-5-5
编辑:Matt J 2019-5-5
nlparci probably wouldn't be applicable if you used lsqnonlin with bound constraints, unless the bounds were inactive at the solution (i.e., you get the same solution as when you don't apply bounds).

类别

Help CenterFile Exchange 中查找有关 Interpolation 的更多信息

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by