I am working with a set of bivariate data arranged into columns labelled 'x' and 'y'. I also have measurements for the error variances corresponding to each observation, labelled 'sx' and 'sy'. An example of the layout is given below:
x sx y sy
1 0.1 2 0.05
Assuming that the measurements for (x,y) correspond to the bivariate normal distribution
$$(\mathbf{x},\mathbf{y})^T \sim N((\mu_1,\mu_2)^T,\Sigma)$$
$$\Sigma=\begin{bmatrix} \tau_{11}+s_{x,i} & \tau_{12} \\ \tau_{12} & \tau_{22}+s_{y,i} \end{bmatrix}$$
where $s_{x,i}$ corresponds to the known error variance of the $i$th $x$ observation etc, how does one go about finding maximum likelihood estimates for the model parameters $(\mu_1, \mu_2, \tau_{11}, \tau_{12}, \tau_{22})$ in R in this case?
Can someone provide code in R which achieves this? The inclusion of the known errors in the covariance matrix is certainly complicating things.