Combine overlapping Gaussian distributions Given the mean and covariance matrix of the pair-wise (multivariate normal) distributions
\begin{align}
P(X_1,X_2) &\sim \mathcal{N}(\mu_{12},\Sigma_{12}),\\
P(X_2,X_3) &\sim \mathcal{N}(\mu_{23},\Sigma_{23}) \quad \text{and}\\
P(X_1,X_3) &\sim \mathcal{N}(\mu_{13},\Sigma_{13})
\end{align}
What is the best estimate for $\mu^\star,\Sigma^\star$ in
\begin{align}
P(X_1,X_2,X_3) &\sim \mathcal{N}(\mu^\star,\Sigma^\star),
\end{align}
the joint distribution of $X_1,X_2,X_3$?

Edit:
Note, that $\mu_{12},\mu_{13},\mu_{23} \in \mathbb{R}^2$ and $\mu^\star \in \mathbb{R}^3$ (similar for the covariance matrices).
If I would have only two one dimensional distributions, $P(X_1) \sim \mathcal{N}(\mu_{1},\Sigma_{1})$ and $P(X_2) \sim \mathcal{N}(\mu_{2},\Sigma_{2})$ I would combine them by concatenating the mean vectors and using $diag(\Sigma_{1}, \Sigma_{2})$ as a diagonal covariance matrix.
 A: This question can have multiple answers depending on who is doing the answering, or, depending on the interpretation placed on it,
it might be unanswerable.
As I understand the question, we are given that $X_1, X_2$, and $X_3$ are pairwise jointly normal random variables with specified mean vectors
and covariance matrices. We are asked for a trivariate joint normal
distribution that is "most consistent" with this information.
As an engineer (and wannabe probabilist), the first thing I would do 
is to check the given information for consistency.  


*

*It must be the case that $\mu_{ij} = (\mu_i, \mu_j)$
where $\mu_i = E[X_i]$, $i = 1,2,3$. So, if the problem claims that $\mu_{12} = (4,10)$ while $\mu_{23} = (5,-3)$, the given information
is inconsistent.

*Similarly, the covariance matrices must be consistent: in 
$\Sigma_{ij}$, the diagonal terms are 
$\sigma_i^2 = \operatorname{var}(X_i)$ and
$\sigma_j^2 = \operatorname{var}(X_j)$. Each variance occurs
in two of the three covariance matrices and it must have the same
value in both occurencse. It must also be true that the
off-diagonal terms in $\Sigma_{ij}$ must be equal, and that
$\Sigma_{ij}$ must have determinant $\geq 0$.

If the above consistency checks are satisfied, then, regardless
  of whether the $X_i$ are pairwise jointly normal or not, the
  mean vector of $(X_1,X_2,X_3)$ is 
  $$\mu = (\mu_1,\mu_2,\mu_3)\tag{1}$$ and
  the covariance matrix is
  $$\Sigma = \left[\begin{matrix}
\sigma_1^2& \sigma_{12}& \sigma_{13}\\
\sigma_{21} & \sigma_2^2 & \sigma_{23}\\
\sigma_{31} & \sigma_{32} & \sigma_3^2\end{matrix}\right]\tag{2}$$
  in which the upper left $2\times 2$ submatrix is the given
  $\Sigma_{12}$, the lower right $2\times 2$ submatrix is the given
  $\Sigma_{23}$, and the "four-corners" submatrix is the given $\Sigma_{13}$.


At this point, the probabilist's answer to the question

Which trivariate joint normal distribution is most consistent
  with the given information?

would be that $\mathcal N(\mu,\Sigma)$ is the only
trivariate normal distribution that is exactly
consistent with the
given information (that the marginal bivariate distributions
are the specified bivariate joint normal distributions). No other trivariate
normal distribution can be as good a match with the given 
information. 
On the other hand, the question might be interpreted as

Which trivariate joint normal distribution is the "best
  approximation" to the actual unspecified joint distribution
  of $(X_1,X_2,X_3)$?

and to my lowly engineering mindset, this question is unanswerable
since neither the meaning of "best" nor the target to be
matched is known;
dyed-in-the-wool probabilists might choose to give a different 
response.  Other readers who disagree with my
interpretation might want to try their hand at
finding the trivariate joint normal distribution that best approximates
$$f(x_1,x_2,x_3) = \begin{cases}
2\phi(x_1)\phi(x_2)\phi(x_3), 
& \text{for}~ x_1 \geq 0, x_2 \geq 0, x_3 \geq 0,\\
&\text{ or}~ x_1<0, x_2< 0, x_3 \geq 0,\\
&\text{ or}~ x_1<0, x_2\geq 0, x_3 < 0,\\
&\text{ or}~ x_1\geq 0, x_2< 0, x_3 < 0,\\
0, &\text{otherwise},\end{cases}$$
where $\phi(\cdot)$ is the standard normal density function.
Note that the $X_i$ are pairwise independent standard normal 
random variables.


On the other hand, a good statistician might take it that the OP
has masses of data (maybe even "Big Data") from which
the OP has drawn the inference that it is reasonable to
assume that the $X_i$ are pairwise jointly normal
random variables, and that $\mu_{12}$, $\mu_{23}$, $\mu_{13}$, $\Sigma_{12}$
$\Sigma_{23}$, $\Sigma_{13}$ are sample means, variances, and covariances
 gleaned from three different data sets
of the form $\{(x_{1,i}, x_{2,i})\colon i = 1,2,\ldots n_1\}$,
$\{(x_{2,j},x_{3,j})\colon j = 1,2,\ldots, n_2\}$ and
$\{(x_{1,k},x_{3,k})\colon k = 1,2,\ldots n_3\}$. Then it
is almost certain that the first entry in $\mu_{12}$
(which is $\frac{1}{n_1}\sum_i x_{1,i}$) is not the same
as the first entry in $\mu_{13}$ (which is $\frac{1}{n_3}\sum_k x_{1,k}$).
Similarly for the other means, variances and covariances.  So,
should the data be merged to get better estimates of what the
means, variances, etc are?  Is there data available of the
form $\{(x_{1,m}, x_{2,m}, x_{3,m})\colon m = 1,2,\ldots N\}$
that might allow some form of estimation as to what the 
 actual joint distribution of $(X_1,X_2,X_3)$ is?
Considerations of this type
lead to the kinds of answers given by Tim and Karel Macek.
A: The problem can be solved as an optimization one. Given $\mu$ and $\Sigma$ for $X_1,X_2,X_3$, you have corresponding marginalizations $\tilde{\mu}_{1,2}$, $\tilde{\Sigma}_{1,2}\dots$. The marginalization of multivariate Gaussian is simple.
To quantify how those match to $\tilde{\mu}_{1,2}$, $\tilde{\Sigma}_{1,2}\dots$, you can use Kulback-Leibler divergence for Gaussians for each of those marginalizations.
Then you have to solve 
$$[\mu^*,\Sigma^*]=\arg\min_{\mu,\Sigma}(KLD_{1,2}+KLD_{2,3}+KLD_{1,3})$$.
A: You are referring to multivariate normal distribution parametrized by a vector of means $\boldsymbol{\mu}$ and covariance matrix $\boldsymbol{\Sigma}$. The MLE estimators for those parameters are sample means for $\mu_i$'s and sample covariance for $\boldsymbol{\Sigma}$.
