nonlinear curve fitting - propagation of data uncertainties to nlpredci confidence intervals?

Hi! I'm trying to fit a nonlinear function that has several parameters, some partially correlated. For simplicity, let's say the function is something like y = x * alpha^x, and the parameter to be fitted is alpha. The data are vectors X and Y. The data have measurement uncertainties that are Gaussian distributed, not correlated, and can have different standard deviations for each x_i or y_i.

I used LSQCURVEFIT to estimate 'alpha' for a given dataset, without applying weights (which ideally I'd like to do, but in R2010a the similar function NLINFIT does not support yet the 'Weights' parameter).

In any case, LSQCURVEFIT finds the best-fit 'alpha' and returns the residuals and a Jacobian. NLPARCI then can compute a 68% confidence interval for 'alpha'.

However, this 68% CI for 'alpha' is just for the unweighted fit itself.

How to get original uncertainties in both X and Y propagated into the 68% CI estimate of 'alpha'?

Answer 1

[Matt J]

How to get original uncertainties in both X and Y propagated into the 68% CI estimate of 'alpha'?

Since alpha, the residuals, and the Jacobian are derived from X and Y, whatever uncertainty was originally in X and Y should already have had its impact on them.

[Jan Wagner]

I'm trying to see that, but am not really convinced yet. My original data are X_i ~ Norm(mu_i,sigma_i) and Y_i ~ Norm(mu_i,sigma_i), where sigma_i!=sigma_j. Perhaps if they were i.i.d. with all sigma_i=sigma_j then NLPARCI might produce the correct parameter 68% CI's.

But if the original distributions are different, shouldn't that either narrow or widen the 68% CI's of the best fit parameters determined by an (unweighted) lsqcurvefit fit?

Btw sorry about the mistake in the title of this post (nlpredci should've been nlparci).

[Matt J]

I doubt the random variation in X_i is accounted for at all. For that, you would have to be doing some sort of Total Least Sqaures approach.

Furthermore, I'm guessing that the CI calculation is based on the following approximation to the Cov matrix of the parameters

 pinv(J)*Cov(Y)*pinv(J).'

where J is the Jacobian and Cov(Y) is the covariance matrix of Y. The question is how the residuals are used by NLPARCI to approximate Cov(Y). Possibly residuals(i) are averaged over i leading to the approximation,

 Cov(Y)=mean(residuals)*eye(N);

However, if you have a better approximation of Cov(Y), you could probably do your own calculation with that instead.