This question is the theoretical counterpart to another question posted on StackOverflow, where I asked about the implementation of the fitting algorithm using Scipy or lmfit libraries for Python. There, it was suggested that I should also ask the question here to have a grasp of the underlying mathematics and perhaps write the code by myself without any libraries.
My question is: How should I proceed if I want to fit a general non-linear model with multiple inputs, multiple outputs and multiple parameters to a dataset that includes covariance matrices for each of the data points? What would be your approach to this? (this can be answered with a least-squares implementation or similar procedures).
To give a more concrete and simple example, suppose I have a 2D vector field, $f:\mathbb{R}^2\longmapsto \mathbb{R}^2$, defined by the function
$$f(x,y;\theta) = {u(x,y) \brack v(x,y)} = {ax+by \brack axy-c}$$
, where $\theta = ( a,b,c )$ is the set of parameters that ought to be adjusted during the fitting procedure, $(x,y)$ are the coordinate positions and $(u,v)$ the components of the vectors.
For some specific parameter values, for example $( a,b,c ) = (-1,7,11)$, the vector field, $f$, looks like this:
Now, suppose I have gathered some data on this vector field, $(x_i, y_i, u_i, v_i)$ for $i = 1,...,N$. Not only that, but my data has some associated uncertainties $(\sigma_{x_i}, \sigma_{y_i}, \sigma_{u_i}, \sigma_{v_i})$ and correlations $(\rho_{xy_i}, \rho_{xu_i}, \rho_{xv_i}, \rho_{yu_i}, \rho_{yv_i}, \rho_{uv_i})$. Thus, for the $i$th datapoint, $(x_i, y_i, u_i, v_i)$, there is an associated covariance matrix
$$ C_i = \begin{pmatrix} \sigma_{x_i}^2 & \sigma_{x_i}\sigma_{y_i}\rho_{xy_i} & \sigma_{x_i}\sigma_{u_i}\rho_{xu_i} & \sigma_{x_i}\sigma_{v_i}\rho_{xv_i}\\\ \sigma_{x_i}\sigma_{y_i}\rho_{xy_i} & \sigma_{y_i}^2 & \sigma_{y_i}\sigma_{u_i}\rho_{yu_i} & \sigma_{y_i}\sigma_{v_i}\rho_{yv_i}\\\ \sigma_{x_i}\sigma_{u_i}\rho_{xu_i} & \sigma_{y_i}\sigma_{u_i}\rho_{yu_i} & \sigma_{u_i}^2 & \sigma_{u_i}\sigma_{v_i}\rho_{uv_i}\\\ \sigma_{x_i}\sigma_{v_i}\rho_{xv_i} & \sigma_{y_i}\sigma_{v_i}\rho_{yv_i} & \sigma_{u_i}\sigma_{v_i}\rho_{uv_i} & \sigma_{v_i}^2 \end{pmatrix} $$
My goal is to have the best fit parameters, $\hat{\theta} = (a_{best}, b_{best}, c_{best})$, for the model, using my data, and considering the full covariance matrix $C_i$ of that data. I also want the covariance matrix of the parameters, because I want to retrieve their uncertainties $(\sigma_{a}, \sigma_{b}, \sigma_{c})$ and correlations $(\rho_{ab}, \rho_{ac}, \rho_{bc})$.
I know only the basics of statistics, so be gentle with me, please.
