This is a reference request. I would like to know how to determine fit parameters along with their uncertainties, if the data have uncertainties in both arguments and measurements.
To be more specific, suppose we have a data set $\big\{\left(x_i,y_i\right)\big\}_{i=1}^n$ along with their uncertainties $\big\{\left(\delta x_i,\delta y_i\right)\big\}_{i=1}^n$. The real numbers $x_i$'s are called the arguments, whereas the real numbers $y_i$'s are their corresponding measurements. The uncertainties in the arguments $\delta x_i$'s and the uncertainties in the measurements $\delta y_i$'s are positive real numbers. The task is to determine the best interpolation $$y=f\left(x;p_1,p_2,\ldots,p_k\right)\,,$$ where $f$ is a given function with parameters $p_1,p_2,\ldots,p_k$ to be determined. The fit function is assumed to be sufficiently nice (e.g., $\mathcal{C}^1$, piecewise smooth, or invertible). Where do I find a method of finding the best fit parameters $\hat{p}_i$'s and their uncertainties $\delta p_i$'s (preferably along with an explanation of why the method should result in a good fit)?
The method I know---the maximum likelihood method---does not take the errors in the arguments $\delta x_i$'s into account. Downloadable codes are also very welcome. If there is a better place to ask this question, please let me know.
Here is a brief description of the maximum likelihood method. Write $\mathbf{p}$ for $\left(p_1,p_2,\ldots,p_k\right)$. The likelihood function of the $i$-th measurement is $$L_i\left(\mathbf{p}\right)=\frac{1}{\sqrt{2\pi\,\delta y_i^2}}\,\exp\Biggl(-\frac{1}{2}\left(\frac{y_i-f\left(x_i;\mathbf{p}\right)}{\delta y_i}\right)^2\Biggr)\,.$$ Then, the total likelihood is then $$L\left(\mathbf{p}\right)=\prod_{i=1}^n\,L_i\left(\mathbf{p}\right)\,.$$ Hence, the maximum likelihood is achieved via implementing the least-squared method on the chi-square value $$\chi^2\left(\mathbf{p}\right):=\sum_{i=1}^n\,\left(\frac{y_i-f\left(x_i;\mathbf{p}\right)}{\delta y_i}\right)^2\,.\tag{*}$$ The uncertainty computation for each parameter $p_j$ is done via finding $$\delta p_j=\sqrt{\sum_{i=1}^n\,\left(\left.\frac{\partial p_j}{\partial x}\right|_{\mathbf{p}=\hat{\mathbf{p}},x=x_i,y=y_i}\,\delta x_i\right)^2+\sum_{i=1}^n\,\left(\left.\frac{\partial p_j}{\partial y}\right|_{\mathbf{p}=\hat{\mathbf{p}},x=x_i,y=y_i}\,\delta y_i\right)^2}\,,\tag{#}$$ where $\hat{\mathbf{p}}$ is a least-square solution to (*).
Surely, (#) does encapsulate the errors in the arguments as well as in the measurements, but the best fit parameters are not sensitive to the uncertainties in the arguments (which I think they should). Maybe if I introduce the likelihood functions of the arguments, namely terms of the form $$\frac{1}{\sqrt{2\pi\,\delta x_i^2}}\,\exp\Biggl(-\frac{1}{2}\left(\frac{x_i-f^{-1}\left(y_i;\mathbf{p}\right)}{\delta x_i}\right)^2\Biggr)\,,$$ provided that $f(\_;\mathbf{p})$ is invertible or at least locally invertible, into the standard likelihood function $L(\mathbf{p})$ above, and employ the same least-square algorithm, then it may be the answer I am seeking. However, I am not certain if this is how data analysts and statisticians work. I understand that the involvement of the inverse function in my approach makes this method very nasty. Hence, I am requesting a reference, in hope that there may be a better method.
May be, you could consider orthogonal distance regression or total least square regression.
The source code of programs ODRPACK and ODRPACK95 are available (for free) from Netlib. I started using the original one almost twenty years ago and it really fitted all my needs. The latest version includes the possibility of bound constraints.
For the method and use of the code, have a look here.