# 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}$$

• Your assumptions about the missing covariances might be mathematically impossible. – whuber Oct 14 '18 at 18:58
• is it not sufficient that the covariance matrix elements satisfy $cov(x_i,x_i)cov(x_j,x_j)>cov(x_i,x_j)^2$? or is there some other way the assumption of zero for the unobserved covariances might be problematic? – zamazalotta Oct 15 '18 at 1:41
• There seem to be other assumptions in what make the "neural stitching" I referred you to, possible. I'm a bit busy for a few days, but your question made me rethink my understanding of this topic, and I'll look into the literature a bit more. Great question, though. – ken Oct 15 '18 at 14:17
• No, that is not sufficient. The determinants of all principal minors must be non-negative. For instance, $$\pmatrix{1&-1/2&0&0\\-1/2&1&-1/2&-1/2\\0&-1/2&1&-1/2\\0&-1/2&-1/2&1}$$ is not a covariance matrix even though it satisfies every one of your pairwise requirements. Notice that its determinant is negative. Indeed, if you replace every $-1/2$ by a constant value that is between $-1$ and $-0.460811$ or between $0.675131$ and $1,$ this matrix will not be a covariance matrix. – whuber Oct 15 '18 at 14:58
• thank you for this counter example. then i will look into possible choices of the unobserved covariance elements that maximize the entropy of the resulting Gaussian – zamazalotta Oct 15 '18 at 15:23

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


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.

• thank you for the detailed answer, I will create a simple example (from a known multivariate normal distribution) and check the likelihood I obtain with the formulation you propose. btw in what literature are such topics covered? does this follow naturally from the definition of the posterior or is there more to it? – zamazalotta Oct 14 '18 at 16:58
• Basically what I used is just the Bayes theorem plus the definition of the independence of datasets, so I'd say any Bayesian statistics textbook should do. – sega_sai Oct 14 '18 at 17:15
• Now that I try to implement this formulation it looks as if you assume $\Sigma_1$ and $\Sigma_2$ have the same dimension (4x4), which is not the case. They each have missing dimensions, one is missing $x_3,x_4$ and the other is missing $x_1$, hence the title of the question. – zamazalotta Oct 15 '18 at 13:05
• You are incorrect. I do not assume that $\Sigma_1$ and $\Sigma_2$ are of the same dimension. If they would be the whole last equation wouldn't make sense because the projection matrices $P_1$ and $P_2$ are not square but rectangular with sizes of 4x2 and 4x3 respectively (4 columns in each). – sega_sai Oct 15 '18 at 15:28
• oops, my apologies... i will try implementing it again – zamazalotta Oct 15 '18 at 15:46

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

• thank you for your interest in my question. I checked the paper and also another review of it and it looks like it wouldn't work in this case as it "requires that the full complexity of the entire matrix is represented in each overlap" [1] Anyhow, given that we never observe cov(x_1,x_3) how can any algorithm deduce (invert) it? To me it looks as if it can only be assumed, and then the only non arbitrary assumption seems to be 0. [1] pillowlab.wordpress.com/2017/10/05/… – zamazalotta Oct 13 '18 at 20:55
• You are right, it would have to be in a covariance matrix which is not full rank. I'm trying to wrap my head around what a covariance matrix that doesn't have full rank looks like. The trivial case with 0s in any element that's connected with say $i$th entry is easy, but I'm trying to picture what sort of arrangement of experiment and observations would yield a low-rank covariance matrix whose diagonal has no 0 in it. – ken Oct 14 '18 at 12:17
• The author of the paper I cited also has a PhD thesis where he applies covariance completion to analysis of neural data where only a subset of neurons are observed at once. reports-archive.adm.cs.cmu.edu/anon/ml2015/CMU-ML-15-106.pdf pg. 41. I can't imagine such data being low rank, so it might be worth taking a look there. – ken Oct 14 '18 at 12:21

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:

1. The covariance matrix is filled with the covariance/variances observed in the experiments
2. If there are repeating elements then they are averaged
3. 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.
4. The optimization ensures that the correlation matrix is symmetric (by design) and positive semidefinite (by ensuring all eigenvalues are positive)
5. 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