Determination of the confidence interval for fitted curves
20 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Torsten Klement
el 1 de Oct. de 2019
Comentada: Adam Danz
el 24 de Oct. de 2019
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
5 comentarios
Respuesta aceptada
Adam Danz
el 22 de Oct. de 2019
Editada: Adam Danz
el 22 de Oct. de 2019
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)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/244002/image.png)
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)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/244003/image.png)
3 comentarios
Más respuestas (1)
Torsten Klement
el 24 de Oct. de 2019
1 comentario
Adam Danz
el 24 de Oct. de 2019
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).
Ver también
Categorías
Más información sobre Solver Outputs and Iterative Display en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!