Solve for Maximum Likelihood Estimate This is the original problem: Suppose that the two-dimensional vectors ($X_1$, $Y_1$), ($X_2$, $Y_2$), ..., ($X_n$, $Y_n$) form a random sample from a bivariate normal distribution for which the means of $X$ and $Y$ are unknown but the variances of $X$ and $Y$ and the correlation between $X$ and $Y$ are known. Find the M.L.E.’s of the means.
(I skipped some steps to get this equation.) Below is the partial derivative equation that I am working on. 
$$\frac{\partial L(\mu_1,\mu_2 )}{\partial \mu_1} = \left [ \frac{1}{(1-\rho ^2)\sigma_1^2}\sum_{i=1}^{n}(x_i-\mu_1)\right ]- \left [ \frac{\rho}{(1-\rho ^2)\sigma_1\sigma_2}\sum_{i=1}^{n}(y_i-\mu_2)\right ] = 0$$
In the equation, $\mu_1$, $\mu_2$ are unknown means for two populations, $\sigma_1$, $\sigma_2$ are known variances for two populations, and $\rho$ is the correlation for bivariate variable $(X, Y)$. $L$ is the likelihood function for $(\mu_1, \mu_2)$.
This question is from a statistics textbook. I want to solve for $\mu_1$ in order to get the MLE (Maximum Likelihood Estimate) of $\mu_1$, but I have trouble in getting rid of $\mu_2$. 
 A: Let us move to a matrix representaion. Define ${\bf{R}}=\begin{pmatrix}
\sigma_1^2 & \rho \sigma_1 \sigma_2\\ 
\rho \sigma_1 \sigma_2 & \sigma_2^2
\end{pmatrix}$ to be the covariance matrix, further define $\underline{\mu}=(\mu_1, \mu_2)^T$, $\underline{x}_i=(x_i, y_i)^T$ and $\underline{x}=(\underline{x}_1^T, \underline{x}_2^T, \dots, \underline{x}_n^T)^T$, so we can write
\begin{equation}
f\left ( \underline{x} \right ) = \frac{1}{\left ( 2 \pi \right )^n \det\left ( R \right )^n} exp \left (  -\frac{1}{2} \sum_{i=1}^{n} \left ( \underline{x}_i-\underline{\mu} \right )^T {\bf{R}}^{-1}\left ( \underline{x}_i-\underline{\mu} \right ) \right ).
\end{equation}
Taking the natural logarithm denoted as $\log$ will yield
\begin{equation}
\log f\left ( \underline{x} \right ) = -nlog\left ( 2\pi \right ) -n \log \left [ \det \left ( R \right ) \right ] - \frac{1}{2} \sum_{i=1}^{n} \left ( \underline{x}_i-\underline{\mu} \right )^T {\bf{R}}^{-1}\left ( \underline{x}_i-\underline{\mu} \right ).
\end{equation}
Deriving with respect to the $2 \times 1$ vector of means, $\underline{\mu}$, will give
\begin{equation}
\frac{\partial }{\partial \underline{\mu}}\log f\left ( \underline{x} \right ) = -\frac{1}{2}\sum_{i=1}^{n}\left [ -2 \underline{x}_i^T {\bf{R}}^{-1} + \underline{\mu}^T {\bf{R}}^{-1} \right ] = 0,
\end{equation}
which gives the result of the well known sample mean:
\begin{equation}
\underline{\hat{\mu}}_{ML} = \frac{1}{n}\sum_{i=1}^{n} \underline{x}_i.
\end{equation}
