Joint probability of multivariate normal distributions with missing dimensions Suppose I conduct two experiments, each measuring a subset of possible parameters. From experiment #1 I measure two parameters and estimate the multivariate normal distribution
$$
\mathbf{X}_1=\left [ x_1,x_2 \right ]
$$
$$
\mathbf{X}_1\ \sim\ \mathcal{N}_1(\boldsymbol\mu_1,\, \boldsymbol\Sigma_1)
$$
$$
\mu_1=[\mu_1^1,\mu_2^1]
$$
$$
\Sigma_1 = \begin{bmatrix}
var(x_1^1) &  \\ 
cov(x_2^1,x_1^1) & var(x_2^1)\\
\end{bmatrix}
$$
In experiment #2 I measure three parameters and build a second multivariate normal distribution
$$
\mathbf{X}_2=\left [ x_2,x_3,x_4 \right ]
$$
$$
\mathbf{X}_2\ \sim\ \mathcal{N}_2(\boldsymbol\mu_2,\, \boldsymbol\Sigma_2)
$$
$$
\mu_2=[\mu_2^2,\mu_3^2,\mu_4^2]
$$
$$
\Sigma_2 = \begin{bmatrix}
 var(x_2^2)& & \\ 
  cov(x_3^2,x_2^2) & var(x_3^2) & \\ 
  cov(x_4^2,x_2^2) & cov(x_4^2,x_3^2) & var(x_4^2) 
\end{bmatrix}
$$


*

*My question is how do I calculate the joint probability distribution that describes the complete space $ \mathbf{X}=\left [ x_1, x_2,x_3,x_4 \right ]$?

*My goal is to use this joint probability distribution to compute the likelihood of a validation set and make model selection.


The formulas for calculating the product of two multivariate pdfs consider the same dimensions, that's why I am confused. 
EDIT: 
I have been thinking about this and here is where I am at:
As Ken puts in his answer, we have not observed $x_1$ and $x_3$ together so we have no estimate for $cov(x_1,x_3)$. So in the absence of this information it looks to me as my best option is to assume $cov(x_1,x_3)=0$ ? 
If this assumption makes sense then can I use the following means and covariances to estimate the product? Notice that I am am "completing" the covariance matrix of experiment #1 with the covariances observed in experiment #2 and vice versa, where the $x_i^j$ denotes the $i$th parameter observed in experiment $j$
$$
\mu_1=[\mu_1^1,\mu_2^1,\mu_3^2,\mu_4^2]
$$
$$
\Sigma_1 = \begin{bmatrix}
var(x_1^1) &  &  & \\ 
cov(x_2^1,x_1^1) & var(x_2^1)& & \\ 
0 &  cov(x_3^2,x_2^2) & var(x_3^2) & \\ 
0 &  cov(x_4^2,x_2^2) & cov(x_4^2,x_3^2) & var(x_4^2) 
\end{bmatrix}
$$
And for experiment #2
$$
\mu_2=[\mu_1^1,\mu_2^2,\mu_3^2,\mu_4^2]
$$
$$
\Sigma_2 = \begin{bmatrix}
var(x_1^1) &  &  & \\ 
cov(x_2^1,x_1^1) & var(x_2^2)& & \\ 
0 &  cov(x_3^2,x_2^2) & var(x_3^2) & \\ 
0 &  cov(x_4^2,x_2^2) & cov(x_4^2,x_3^2) & var(x_4^2) 
\end{bmatrix}
$$
 A: Assuming that the one can interpret the results of the experiments as  that the likelihood of the data is 
$P(D_1|{\bf x}) \propto N([x_1,x_2]|\mu_1,\Sigma_1)$
$P(D_2|{\bf x}) \propto N([x_2,x_3,x_4]|\mu_2,\Sigma_2)$
Then one can write down the posterior for two experiments combined (assuming data are independent)
$P({\bf x}|D) \propto P(D_1|{\bf x}) P (D_2|{\bf x}) \pi({\bf x})$
To determine the exact distribution here is easy, if we define two projection matrices  $P_1$ and $P_2$, such that 
$[x_1,x_2] = P_1 \bf{x}$ and
$[x_2,x_3,x_4] = P_2 \bf{x}$
Then the log posterior is 
$ -2 \log(P) = (P_1 {\bf x} - \mu_1)^T \Sigma_1^{-1} (P_1 {\bf x} - \mu_1) + (P_2 {\bf x} -\mu_2)^T \Sigma_2^{-1}(P_2 {\bf x} - \mu_2)$
Rearranging we get 
$$ {\bf x}^T (P_1^T \Sigma_1^{-1} P_1 + P_2^T \Sigma_2^{-1} P_2) {\bf x} - 2 (\mu_1^T \Sigma_1^{-1} P_1 + \mu_2^T \Sigma_2^{-1} P_2 ) {\bf x} + \mu_1^T \Sigma_1^{-1} \mu_1 + \mu_2^T \Sigma_2^{-1} \mu_2$$
Therefore the posterior on $\bf{x}$ is 
$$\mathcal{N}((P_1^T \Sigma_1^{-1} P_1 + P_2^T \Sigma_2^{-1}P_2)^{-1} (P_1^T\Sigma_1^{-1}\mu_1+P_2^T\Sigma_2^{-1}\mu_2), (P_1^T \Sigma_1^{-1} P_1 + P_2^T \Sigma_2^{-1} P_2)^{-1})$$
EDIT: 
Apparently, I misunderstood what the original problem was, so I was solving something different. For the actual problem in hand, as ken pointed you don't seem to have enough information, however, you can get the results if you assume some sort of prior on the covariance matrix (i.e. Wishart). Then 
You can write the likelihoods of two datasets (where $P_1$, $P_2$ are projection matrices) 
$$P(D_1|\mu,\Sigma) = N(D_1|P_1 \mu,P_1\Sigma P_1^T)$$
$$P(D_2|\mu,\Sigma) = N(D_2|P_1 \mu,P_2\Sigma P_2^T)$$
 Then you can combine those with the Wishart prior and maximize.
$$P(\mu,\Sigma|D) \propto \pi(\mu)\pi(\Sigma)P(D_1|\mu,\Sigma) P(D_2|\mu,\Sigma)$$
 I am sure there is an analytical expression that comes out of this,  but I don't have time to derive it. 
A: Assuming $cov(x_1,x_3)=0$ is unjustified and leads to inconsistencies in your co-variance matrix.
Instead you can use MMSE estimators calculated based on the available data to fill in unknown elements. Covariance between $x_1$ and $x_2$ can be calculated from the first measurement and between $x_2$ and $x_{3,4}$ from the second measurement.
I suspect that this would be the result of analysis by @sega_sai if it is completed.
I assume a zero mean for variables when I write MMSE estimators. This however does not change the result of the co-variance.
MMSE estimator of $X_1$ given $X_2$ can be calculated as:
$$
\hat x_1=\frac{cov(x_1,x_2)}{var(x_2)} x_2
$$
$$
\hat{Cov}(x_1,x_3)=Cov(\hat x_1, x_3)=\frac{cov(x_1,x_2) cov(x_2,x_3) }{var(x_2)} 
$$
I did not have time to prove if this always generates a positive definite matrix but I tried some random co-variance matrices and it does reproduce positive definite matrices.
clc
for i=1:1000
A=randn(3,3);
C=A'*A;
C_=C;

c_=C(1,2)*C(2,3)/C(2,2);
C_(3,1)=c_;
C_(1,3)=c_;
if det(C_)<0
    disp('Negative det');
    det(C)
    C
    C_
    i
    break
end
end

for i=1:10000
A=randn(4,4);
C=A'*A;
C_=C;


c_=C(1,2)*C(2,3)/C(2,2);
C_(3,1)=c_;
C_(1,3)=c_;

c_=C(1,2)*C(2,4)/C(2,2);
C_(4,1)=c_;
C_(1,4)=c_;
if det(C_)<0
    disp('Negative det');
    det(C)
    C
    C_
    i
    break
end
end
i

A: Really, you want to recover the $4 \times 4$ covariance matrix that you would've gotten had you measured all $[x_1, x_2, x_3, x_4]$ jointly.  But you didn't.  You measured a subset of variables that happened to overlap.  (overlap is very important for this problem.  Had you measured $[x_1, x_2]$ and $[x_3, x_4]$ - then forget it, no way to recover $4 \times 4$ covariance.)  The idea is that you have matrix components X that were measured from the 1st experiment, and matrix components Y that are measured from the 2nd component.   
$
\left(
\begin{matrix}
X & X & \_ & \_\\
X & (XY) & Y & Y \\
\_ & Y & Y  & Y \\
\_ & Y & Y & Y\\
\end{matrix}
\right)
$
What you're asking is, is there a way to fill in the entries $\_$ in the above covariance matrix, and it is great question, and also not trivial at all.
I don't know all the details, but there is a paper that deals with this problem.  It may be of help for you.
http://papers.nips.cc/paper/5467-deterministic-symmetric-positive-semidefinite-matrix-completion.pdf
A: I would like to thank everyone who contributed to this discussion, especially @sega_sai, @ken for taking the time to implement and provide critical review.
As a result of this fruitful discussion I have converged to the following solution:


*

*The covariance matrix is filled with the covariance/variances observed in the experiments

*If there are repeating elements then they are averaged

*The missing elements are estimated by constrained maximization of the log determinant of the covariance matrix (i.e maximizing the entropy). For this purpose the covariance matrix is converted to a correlation matrix (as all diagonal elements are observed) giving upper and lower bounds of [-1 1] for the unknown off-diagonal elements.

*The optimization ensures that the correlation matrix is symmetric (by design) and positive semidefinite (by ensuring all eigenvalues are positive)

*The resulting correlation matrix is then converted back to covariance and used to estimate the likelihood of an independent dataset. This likelihood is compared to the likelihood of the true (generating) model and the model obtained using Wishart priors (as per @sega_sai).



=== Case 1 ====
True model logL:      -594.89
Wishart prior logL:   -605.75
max Determinant logL: -596.41

=== Case 2 ====
True model logL:      -525.53
Wishart prior logL:   -558.85
max Determinant logL: -531.41

=== Case 3 ====
True model logL:      -513.15
Wishart prior logL:   -962.17
max Determinant logL: -518.87

For the three cases (provided by  @ken) the determinant maximization gives significantly better likelihoods. That is why I am accepting my own answer but awarding the bounty to @sega_sai for his efforts in implementing the alternate solution.
I should also note that @Hooman's MMSE implementation gives the same results for this example, but the formulation requires a common dimension to be observed in the two experiments (in this case $x_2$). The determinant maximization,however, does not require this and therefore is a more general solution.
The Matlab implementation of the determinant maximization method can be accessed here: https://pastebin.com/UTkLUUrv
