Joint distribution of two multivariate normal distributions

Question

If we define 2 independent variables $Y_1$ and $Y_2$ as follows: \begin{align} Y_1 &= (Y_{11},Y_{12},Y_{13})^T \sim\mathcal N_3(\mu_1,\Sigma_{11}), \\ Y_2 &= (Y_{21},Y_{22})^T \sim\mathcal N_2(\mu_2,\Sigma_{22}) \end{align} where, \begin{align} \mu_1 &= (2, 2, 2)^T &\Sigma_{11} &= \left[\begin{array}{ccc} 3 &1 &0 \\ 1 &2 &0 \\ 0 &0 &3 \end{array}\right] \\ \mu_2 &= (3, 4)^T &\Sigma_{22} &= \left[\begin{array}{cc} 4 &2 \\ 2 &4 \end{array}\right] \end{align}

Then how can I find the joint distribution of $Y_{11}-Y_{13}+Y_{22}$ and $Y_{21}-Y_{12}$?

I know its a simple question but I could find if it was asked for $Y_1-Y_2$ or something. How am I supposed to solve it when it is like that?

Answer 1

Assuming that the word independent in the opening statement is used in the way that probabilists use the word and not in the sense of independent versus dependent variable as is common in regression analysis, the joint distribution of the five random variables $Y_{11}, Y_{12}, Y_{13}, Y_{21},Y_{22}$ is the product of the joint distributions of $Y_{11}, Y_{12}, Y_{13}$, and $Y_{21},Y_{22}$, both of which are multivariate normal. This $5$-variate joint distributions is also a multivariate normal distribution in which the mean vector is just the concatenation $(\mu_1, \mu_2)^T$ of the two mean vectors and the covariance matrix is $$\Sigma = \left[\begin{matrix}\Sigma_{11} & 0\\0 & \Sigma_{22}\end{matrix}\right].$$ Thus, the joint distribution of $Y_{11}-Y_{13}+Y_{22}$ and $Y_{21}-Y_{12}$ is a bivariate normal distribution which can be found by the standard methods involving setting up a linear transformation mapping $(Y_{11}, Y_{12}, Y_{13}, Y_{21},Y_{22})$ to $Y_{11}-Y_{13}+Y_{22},Y_{21}-Y_{12})$ and doing matrix calculations. More simply, the means and variances of $Y_{11}-Y_{13}+Y_{22}$ and $Y_{21}-Y_{12}$ as well as their covariance can be computed more directly and used in writing down the mean vector and covariance matrix of this bivariate normal distribution.

Answer 2

Given that $Y_{1}$ and $Y_{2}$ are independent, we have that

$$ \left[\array{Y_{11} \\ Y_{12} \\ Y_{13} \\ Y_{21} \\ Y_{22}}\right] \sim MVN\left(\left[\array{2\\2\\2\\3\\4}\right],\left[\array{3 & 1 & 0 & 0 & 0\\1 & 2 & 0 & 0 &0\\0&0&3&0&0\\0&0&0&4&2\\0&0&0&2&4}\right]\right) $$

Let

$$ \begin{array}{rcl}X_1 & = & Y_{11}-Y_{13}+Y_{22}\\ X_2 & = & Y_{21}-Y_{12}\end{array} $$

As $Y_{11},Y_{12},Y_{13},Y_{21},Y_{22}$ are jointly normal, the linear combinations $Y_{11}-Y_{13}+Y_{22}$ and $Y_{21}-Y_{12}$ are normally distributed. It also follows that as any linear combination of $X_{1}$ and $X_{2}$ is a linear combination of $Y_{11},Y_{12},Y_{13},Y_{21},Y_{22}$ so must $X_{1}$ and $X_{2}$ be jointly normal.

All that remains is to determine the mean and covariance of $X_{1}$ and $X_{2}$. Given the linearity of expectations, the mean is trivial to calculate:

$$ \begin{array}{rcl} E[X_1] &=& E[Y_{11} - Y_{13} + Y_{22}]\\ &=& E[Y_{11}] - E[Y_{13}] + E[Y_{22}]\\ E[X_2] &=& E[Y_{21} - Y_{12}]\\ &=& E[Y_{21}] - E[Y_{12}] \end{array} $$

The covariance is equally straightforward yet tedious:

$$ \begin{array}{rcl} Cov[X_1,X_1] &=& Cov[Y_{11},Y_{11}] + 2 \times Cov[Y_{11},-Y_{13}+Y_{22}] + Cov[-Y_{13}+Y_{22},-Y_{13}+Y_{22}]\\ &=& Cov[Y_{11},Y_{11}] - 2 \times Cov[Y_{11},Y_{13}] + 2 \times Cov[Y_{11},Y_{22}] + Cov[Y_{13},Y_{13}] - 2 \times Cov[Y_{13},Y_{22}] + Cov[Y_{22},Y_{22}]\\\\ Cov[X_2,X_2] &=& Cov[Y_{21},Y_{21}] - 2 \times Cov[Y_{12},Y_{21}] + Cov[Y_{12},Y_{12}]\\\\ Cov[X_1,X_2] &=& Cov[Y_{11},Y_{21}-Y_{12}] + Cov[-Y_{13}+Y_{22},Y_{21}-Y_{12}]\\ &=& Cov[Y_{11},Y_{21}] - Cov[Y_{11},Y_{12}] - Cov[Y_{13},Y_{21}] + Cov[Y_{13},Y_{12}] + Cov[Y_{22},Y_{21}] - Cov[Y_{22},Y_{12}] \end{array} $$

Fortunately many of these terms are zero.

Given the tedious nature of the calculations you can do a simple Monte Carlo simulation to check your answers. Here is some R code for achieving that:

# Include MASS library for mvrnorm for generating multivariate normally distributed samples
library(MASS)

generateSamples <- function(N)
{
  # Generate N samples from Y1 and Y2 with the given mean vectors and covariance matrices
  Y1 <- mvrnorm(mu=rep(2,3),Sigma=matrix(c(3,1,0,1,2,0,0,0,3),nrow=3,ncol=3),n=N)
  Y2 <- mvrnorm(mu=c(3,4),Sigma=matrix(c(4,2,2,4),nrow=2,ncol=2),n=N)

  # Calculate X1 and X2
  X1 <- Y1[,1] - Y1[,3] + Y2[,2]
  X2 <- Y2[,1] - Y1[,2]

  cbind(X1,X2)
}

# Generate 100000 samples from X1 and X2
mySample <- generateSamples(100000)

# Empirical mean vector
mu <- colMeans(mySample)

# Empirical covariance matrix
Sigma <- cov(mySample,mySample)