How to get the confidence intervals of regression coefficients in nlinfit?

I used nlinfit or lsqcurvefit to do non-linear fit.

 [beta,R,J,CovB,MSE,ErrorModelInfo] = nlinfit(X,Y,@myfun,beta0)

How to get the confidence intervals of regression coefficients 'beta'?

I'm sorry it might be just a math problem, but I really want to implement it in matlab code.

Thank you very much!

Answer 1

[Adam Danz]

I looked into this somewhat recently and found two main strategies for computing confidence intervals of the coefficients.

1. using nlparci()
2. computing them directly by calculating the covariance matrix

Both methods use the jacobian and residual outputs and for my data there are very tiny differences in the results between the two methods.

Here are some sources I stored within the footnotes of my code and I recommend browsing through them to determine the best method for your problem.

Approximation of covar matrix of the params

• https://stats.stackexchange.com/questions/231868/relation-between-covariance-matrix-and-jacobian-in-nonlinear-least-squares
• https://www.mathworks.com/matlabcentral/answers/51136-calculate-uncertainty-for-fitted-parameter-from-least-squares-fit

Both covar matrix and nlparci methods

• https://www.mathworks.com/matlabcentral/answers/155477-standard-deviation-of-parameters-in-lsqcurvefit

How covar matrix is used to get CI

• https://stats.stackexchange.com/questions/241449/matrix-and-regression-model?utm_medium=organic&utm_source=google_rich_qa&utm_campaign=google_rich_qa

One method to approximate CIs using linear example, explains how to do with nonlinearly.

• https://www.mathworks.com/matlabcentral/answers/94393-what-is-the-algorithm-used-in-nlparci-in-the-statistics-toolbox

nlparci can be used with lsqnonlin outputs (my project uses lsqnonlin but the outputs are similar for your project)

• https://stackoverflow.com/questions/28212335/how-do-i-calculate-in-matlab-the-95-confidence-interval-with-lsqcurvefit
• http://www.44342.com/matlab-f582-t49885-p1.htm#wrap
• https://www.coursehero.com/file/25475830/MATLAB-TUTORIAL-WEEK-9-jEREMYLEE-44938769docx/
• https://www.mathworks.com/matlabcentral/answers/166726-is-it-statistically-correct-using-the-jacobian-provided-by-lsqnonlin-as-input-for-nlparci-to-get-par

If you go with the nlparci() method, it's as easy as pluggin in your outputs...

   CI = nlparci(betaEstimates,residuals,'jacobian',jacobian);

If you go with the covar matrix method, here a snippet from my code that should be fairly self-explanatory.

    df = length(residuals) - numel(betaEstimates);   %degrees of freedom
    alpha = 0.05;                                    %for 95% conf int
    crit = tinv(1-alpha/2,df);                       %critical value
    covm = inv(jacobian'*jacobian) * var(residuals); %covariance matrix
    covmIdx = sub2ind(size(covm),1:size(covm,1),1:size(covm,2));  %indices of the diag of covm
    CI = nan(numel(betaEstimates),2);

I'm refraining for making a recommendation partially due to what @SK mentioned; there isn't a standard solution to getting CIs for the beta estimates, and also because I'm exploring these methods as a novice (and welcome feedback).

Answer 2

[SK]

There is no fixed answer to this, unless myfun is known and is nice enough to be analysed theoretically. In general, I would do the following:

1. Look at the residuals, Y - myfun(X, beta0). (Here beta0 are your coefficient estimates).

2. Try to determine the distribution of the residuals.

3. Generate many realizations of Y using Y = myfun(X, beta0) + R, where R is generated randomly according to the distribution found in (2). To each realization, do an nlinfit and find the vector beta0_y.

4. Look at the distribution of beta0_y. It may be identifiable as a known distribution or not. Regardless, you can find confidence intervals from it.