Derive multivariate from bivariate normal distribution

Question

Could anyone help me on the following.

Let $K$ and $M$ be integers so that $K\geq3$ and $2\leq M < K$. Let $\boldsymbol{X}=(X_1, ..., X_K)^T$ be a random vector, $\boldsymbol{\mu}$ be a $K\times 1$ vector, and $\Sigma_{K\times K}$ be a positive definite matrix.

We all know that if $\boldsymbol{X}\sim N(\boldsymbol{\mu}, \Sigma)$ then every $M-$variate marginal distribution is also normal.

Suppose inversely that for every $M-$variate marginal distribution is normal, i.e., \begin{equation*} (X_{i1}, ..., X_{iM})^T\sim N(\boldsymbol{\mu^{'}}, \Sigma^{'}) \end{equation*} where $\boldsymbol{\mu^{'}}, \Sigma^{'}$ is obtained by keeping only corresponding rows and columns of $\boldsymbol{\mu}, \Sigma$, respectively. Does this lead to the following \begin{equation*} \boldsymbol{X}\sim N(\boldsymbol{\mu}, \Sigma). \end{equation*}

Thank you so much!

Answer 1

I'm afraid this is not generally true. A counterexample is afforded by emulating a standard example of a non-normal bivariate distribution whose marginals are normal: erase all probability from (say) the second and fourth quadrants, as shown in the upper right example of Cardinal's answer.

Extend this to the case $M=2,K=3$ in a similar way: beginning with a standard trivariate Normal distribution, erase all probability in the octants where an odd number of $X_1,X_2,X_3$ are negative. (Dilip Sarwate provides the details in this answer from several years ago.)

This R example shows how to generate data from this distribution.

n <- 1e4 # Approximately twice the desired sample size
K <- 3   # Must be 3 or greater
Y <- matrix(rnorm(K*n), n, dimnames=list(NULL, paste0("X.", 1:K)))
X <- Y[rowSums(Y >= 0) %% 2 == 0, ]

The last line removes every vector in which an odd number of the components are negative.

A picture should convince you that this works (and will lead easily to a rigorous demonstration that all three $M$-variate marginals are standard Normal). Here is a scatterplot matrix:

pairs(X, pch=16, col="#00000003")

Nevertheless, this is clearly not $K$-variate normal, because the variable determined by counting the positive components does not have a Binomial$(1/2,3)$ distribution:

table(rowSums(X >= 0))

   0    2 
1235 3778 

---

This approach is not special to the case $K=3,M=2:$ it generalizes to all $M$ and $K.$ I want to convince you all the multivariate marginals remain Normal, perhaps by simulating data as above. The challenge in a simulation is to perform an adequate test of multivariate normality in all the margins. One way is to check that random linear combinations are Normal.

For $K=4,$ I conducted that test with 500 random linear combinations for each multivariate marginal, using a chi-squared test based on binning the results into 20 quantiles of the standard Normal distribution. For comparison, I did the same thing with a sample from a multivariate standard Normal distribution of the same size. To compare the results, here are histograms of the chi-squared p-values. The column labels are the columns determining the marginal distributions. The "fake" multivariate Normal is plotted along the top row while the reference ("real") multivariate Normal is plotted along the bottom row.

We expect these histograms to be approximately uniform for a truly $K$-variate Normal distribution. They depart from this a little bit because the linear combinations are not independent of each other. However, by comparing the frequencies of lower p-values in each column but the last, it is abundantly clear that according to these tests the "fake" $K$-variate distribution looks just as Normal for all $M$-variate marginals. Because the frequency of low p-values is so much greater for $K=M$ (rightmost column) it does not look at all like it's $K$-variate Normal, even though all the $M$-variate marginals for $2\le M \lt K$ do appear Normal.

---

For those who would like to experiment, here is the full R code used to generate the examples and figures.

set.seed(17)
n <- 1e4
K <- 3  # Must be 3 or greater
Y <- matrix(rnorm(K*n), n, dimnames=list(NULL, paste0("X.", 1:K)))
X <- Y[rowSums(Y >= 0) %% 2 == 0, ]
Y <- Y[1:nrow(X), ]
Y <- Y * matrix(sample.int(2, length(Y), replace=TRUE)*2 - 3, ncol=K)

pairs(X, pch=16, col="#00000003")

table(rowSums(X >= 0)) # "Fake" multivariate normal distribution
table(rowSums(Y >= 0)) # "Real" MVN distribution
#
# To verify M-variate Normality, study all M-subsets of {1,2,...,K} for
# normality.  A quick check is to take random linear combinations and test
# them for normality.
#
cutpoints <- qnorm(cutpoints.p <- seq(0, 1, by=0.05))
combine <- function(l) if("list" %in% class(l)) do.call("rbind", l) else l

n.tests <- 500
models <- list(`Non-normal`=X, Normal=Y)
df <- combine(lapply(names(models), function(Z.name) {
  Z <- models[[Z.name]]
  combine(lapply(2:K, function(M) {
    combine(apply(combn(K, M), 2, function(j) {
      Margin <- Z[, j]
      p <- sapply(1:n.tests, function(i) {
        y <- rnorm(M)
        z <- Margin %*% (y / sqrt(sum(y^2)))
        chisq.test(table(cut(z, cutpoints)))$p.value
      })
      data.frame(p=p, Model=Z.name, Dimension=M, Margins=paste(j, collapse=","))
    }))
  }))
}))
#
# Plot histograms of the p-values.
#
library(ggplot2)
ggplot(df, aes(p)) +
  geom_histogram(aes(fill=Margins), show.legend=FALSE, binwidth=0.05, boundary=0) +
  facet_grid(Model ~ Margins)

Answer 2

This is an interesting question. I will provide a formal counterexample.

So let us construct a counterexample to show that this does not hold in general. Let us assume $K=3$ and let us consider two independent normal variables $X_1$ and $X_2$ satisfying the standard normal distribution $X_1,X_2\propto \mathcal{N}(0,1)$.

Let us construct the variable $X_3$ as follows

\begin{equation} X_3=\begin{cases} X_2 & for\qquad X_2\geq X_1\\ -X_2 & for\qquad X_2<X_1 \end{cases} \end{equation}

1. Joint distribution of $X_3$ and $X_1$, $P(X_3,X_1)$: joint Normal and independent. Because the realization of $X_2$ does not depend on $X_1$. Reflections at the origin do not change this result.

2. Joint distribution of $X_3$ and $X_2$, $P(X_3,X_2)$: joint Normal and dependent. The threshold construction does not destroy normality as we are integrating out $X_1$ and $X_1$ and $X_2$ are independent by assumption.

3. Joint distribution of $X_1$ and $X_2$: independent and normal by assumption

But clearly, the joint distribution $P(X_1,X_2,X_3)$ can not be expressed as a multivariate Gauss Distribution.

Of course, by assumption! You assume

Suppose inversely that for every M−variate marginal distribution is normal

Hence, the full K-variate marginal distribution must also be normal, isint it? >>>Otherwise it would conflict with your above assumption.