Determination of the confidence interval for fitted curves

Hello,

I have measurement data and I'm fitting it against a differential equation system. In addition, I would like to calculate the two curves that include the measurements within the confidence interval.

In order to do this I have tried the following approach:

[Opts,resnorm,resid,exitflag,output,lambda,J] = lsqcurvefit(@main_fun,Starts,x,y,l_lim,u_lim);
Ci = nlparci(Opts,resid,'jacobian',J);

'Opts' contain two parameters that are to be fitted. In 'Ci' I also get as expected an upper and a lower limit for these values.

The limits are clearly too close to the values only for my sensation.

I have attached a figure of my results and a figure of my expectations.

I actually thought I was on the right track. Did I make a logical mistake?

Thank you and kind regards

Torsten

[Adam Danz]

Hello Torsten , could you show us how you're plotting the prediction intervals produced by nlpredci?

A quick note on nlparci:

Since you're using upper and lower bounds in the fit, in order to compute the confidence intervals of the parameter estimates you'd need to show that the bounds aren't being hit during the fit. If the bounds are being hit, you cannot use nlparci(). That topic and alternative proposals are discussed here and here.

[Torsten Klement]

Hallo Adam,

thanks for the input.

Here are my scripts i tested so far. I wasn't able to solve the problem.

also:

if there is a smoother solution instead of the usage of this mainFunc i would be happy to learn this. But this works fine for all my problems with these kind of problem.

kind regards ,

Torsten

[Adam Danz]

It looks like you're plotting the prediction intervals correctly.

I increased your noise parameter to r = 0.8 in order to see more variation. Then I computed the confidence interval of your parameter estimate (not the prediction interval) in order to compare it to your prediction interval.

% The confidence interval on your parameter estimate is 
CI = nlparci(k_opt,resid,'jacobian',J);
% CI = [0.035 0.052]

% which means the k_opt parameter may vary between [0.035 0.052] within the 95% CI
% Let's plot the function with k_opt parameter varying between those bounds
figure(); 
plot(t_data,Y_data(:,1),'xk') ; 
hold on
param = linspace(CI(1),CI(2),20);
arrayfun(@(x)plot(t_data,mainFunc(x, t_data),'-'),param)
% Show the best-fit line
plot(t_data, mainFunc(k_opt, t_data), 'k--','linewidth',3)

% Now let's overlay the prediction interval
plot(t_data,[lower,upper],'b--', 'linewidth',2)

As you can see, the prediction intervals mostly agree with the range of parameter estimates within the 95% CI.

Additionally, I fit your data using lsqnonlin() and got the same results and same residuals.

[x,resnorm,residual,exitflag,output,lambda,jacobian] = lsqnonlin(@(x)mainFunc(x,t_data)-Y_data, k_start); 

[Torsten Klement]

Hi Adam Danz,

thank you very much for your work on this topic.

But shouldn't 95% of the generated data points be within the limits of the interval?

kind regards

Torsten

[Adam Danz]

I see now. See my answer below.

Answer 1

[Adam Danz]

@Torsten Klement, (continuing from comments under the question)

There are 2 types of prediction intervals.

Confidence interval of the function value at x

In your current prediction interval calculation, you are computing the confidence interval for the estimated curve (function value) at each observation (x). To demonstrate, the 95% confidence interval of your paramter estimate (using nlparci) is CI = [0.035 0.052] which means the k_opt parameter may vary between [0.035 0.052] within the interval. The plot below shows a range of possible curves when the k_opt parameter is varied between those bounds. Then I add the 95% confidence interval for the estimated curve using your line of code (nlpredci). As you can see, it's roughly the bounds of the varied curve. See this example in the documentation which produces a similar result.

% The confidence interval on your parameter estimate is 
CI = nlparci(k_opt,resid,'jacobian',J);
% CI = [0.035 0.052]

% which means the k_opt parameter may vary between [0.035 0.052] within the 95% CI
% Let's plot the function with k_opt parameter varying between those bounds
figure(); 
plot(t_data,Y_data(:,1),'xk') ; 
hold on
param = linspace(CI(1),CI(2),20);
arrayfun(@(x)plot(t_data,mainFunc(x, t_data),'-'),param)
% Show the best-fit line
plot(t_data, mainFunc(k_opt, t_data), 'k--','linewidth',3)

% Now let's overlay the CI of the function
[ypred,delta] = nlpredci(@mainFunc,t_data,x,residual,'Jacobian',jacobian)
lower = ypred-delta;
upper = ypred+delta;
plot(t_data,[lower,upper],'b--', 'linewidth',2)

Prediction for new observations at x

To compute the prediction interval of new observations (that should always be larger than the confidence interval of the data), you need to use the PredOpt-observeration flag.

[ypred,delta] = nlpredci(@mainFunc,t_data,x,residual,'Jacobian',jacobian,'PredOpt','observation')
lower = ypred-delta;
upper = ypred+delta;
plot(t_data,[lower,upper],'r:', 'linewidth',2)

[John D'Errico]

Thank you Adam, I was reading the initial comments on the question itself, and I was going to make exactly your response. Then I saw your answer, so you saved me some time.

[Adam Danz]

I wish we had some kind of indicator that an answer or comment is being written. Glad you caught it before investing time!

[Torsten Klement]

Thanks again. Perfect what i tried to do.

Answer 2

[Torsten Klement]

Actually, this question is answered. As a bonus it would be great if someone could explain how to calculate the delta values from the Jacobi Matrix and the residuals. Unfortunately I can't find a mathematical description for it.

thanks again everybody.

Regards,

Torsten

[Adam Danz]

Are you describing confidence intervals of the parameter estimates? Here's how to compute them using the Jacobian and residuals. Be sure to read both answers there (mine and Matt J's).

https://www.mathworks.com/matlabcentral/answers/460386-can-i-use-the-jacobian-provided-by-lsqnonlin-to-compute-the-confidence-intervals-using-nlparci#answer_373690