Suppose I have to fit a vector field $\vec{y}_i$ with a vector valued function $\vec{f}_i = \vec{f}(\vec{x}_i; \vec{p})$, so both $\vec{y}_i, \vec{f}_i \in \mathbb{R}^m$, where $\vec{x}_i \in \mathbb{R}^n$, and the index $i$ indicates the i-th datapoint. We are interested in the values of the parameters $\vec{p} = (p_1, \ldots, p_k)$; these can be obtained by minimizing the least squares

$\chi^2=\sum_{i=1}^N \Big\vert \frac{\vec{f}(\vec{x}_i) - \vec{y}_i}{\vec{\sigma}_i} \Big\vert^2 = \sum_{i=1}^N \sum_{j=1}^m \Big( \frac{f_j(\vec{x}_i) - y_{ji}}{\sigma_{ji}} \Big)^2$

My question: what is the correct way to estimate the uncertainties in the parameters $\vec{p}$ of this model?

What I have tried: For a scalar function $g_i = g(\vec{x}_i;\vec{p})$, we typically approximate the covariance matrix as

$C = (J_{\vec{p}}^{T} W J_{\vec{p}})^{-1}$,

where $J_{\vec{p}}=\nabla_{\vec{p}} g(\vec{x}_i;\vec{p})$ is the jacobian of $g_i$ with respect to the parameters in the model, and $W$ is a diagonal matrix with the values of $1/\sigma_i^2$ on its diagonal. This way, $J_{\vec{p}}$ is a $k \times N$ matrix, and $W$ is $N \times N$ ($N$ is the number of datapoints).

I therefore assumed the natural way to generalize this to the vector model $\vec{f}_i$, is to make $J$ and $W$ higher dimensional tensors, such that $J_{\vec{p}}$ is $k \times m \times N$, and W is $m\times N \times N \times m$. This way, if I dot product them together as $J_{\vec{p}}^{T} W J_{\vec{p}}$ appropriately, the result will still be a $k \times k$ matrix as it should be, which after inversion should be the covariance matrix. But the results do not quite match my expectation (though they are close) so I want to make sure this is the right way to deal with this. Is my generalization correct, or if not, what should it be?


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.