# Fitting a model with multiple inputs, multiple outputs, multiple parameters, and covariance matrices for each data point

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.