NMST539 | Lab Session 3

Multivariate Normal Distribution

(joint, marginal, and conditional distributions)

Summer Term 2021/2022 | 01/03/22

source file (UTF8 encoding)

Outline for the third NMST539 lab session

multivariate normal distribution (joint, marginal, and conditional;
linear combinations of random vectors with multivariate normal distribution;
random vector properties: normality vs. independence/dependence;

The R-software is available for download from the website: https://www.r-project.org

A user-friendly interface (one of many): RStudio.

Manuals and introduction into R (in Czech or English):

Bína, V., Komárek, A. a Komárková, L.: Jak na jazyk R. (PDF súbor)
Komárek, A.: Základy práce s R. (PDF súbor)
Kulich, M.: Velmi stručný úvod do R. (PDF súbor)
De Vries, A. a Meys, J.: R for Dummies. (ISBN-13: 978-1119055808)

Theoretical/statistical references

Hardle, W. and Simar.L.: Applied Multivariate Statistical Analysis. Springer, 2015
Mardia, K., Kent, J., and Bibby, J.:Multivariate Analysis, Academic Press, 1979.

1. Conditional normal distribution

Let us start by considering a simple two-dimensional normal distribution of a random vector $\Big(\begin{array}{x}X_{1}\\X_{2}\end{array}\Big) \in \mathbb{R}^2$. This is usually formally denoted as

$\Big(\begin{array}{x}X_{1}\\X_{2}\end{array}\Big) \sim N_{2}\left(\boldsymbol{\mu} = \Big(\begin{array}{c} \mu_{1} \\ \mu_{2}\end{array}\Big), \Sigma = \left( \begin{array}{cc} \sigma_{1}^{2} & \sigma_{12} \\\sigma_{21} & \sigma_{2}^{2} \end{array} \right) \right)$,

where $\boldsymbol{\mu} \in \mathbb{R}^2$ is the vector of the mean values and $\Sigma$ is the variance-covariance matrix, which is positive definite and symmetric, thus $\sigma_{12} = \sigma_{21}$. The correspoding density function (of the two dimensional normal distrubution) is given by the expression

$\large{f(\boldsymbol{x}) = \frac{1}{2 \pi |\Sigma|^{1/2}} exp\Big\{ -\frac{1}{2} (\boldsymbol{x} - \boldsymbol{\mu})^{\top} \Sigma^{-1} (\boldsymbol{x} - \boldsymbol{\mu}) \Big\},}$

for an arbitrary element of a two dimensional real plane, thus, $\boldsymbol{x} = (x_{1}, x_{2})^{\top} \in \mathbb{R}^{2}$.

This density can be used to derive the marginal distrubution of the random variables $X_{1}$ and $X_{2}$ or the conditional distrubution of $X_{1}$ given $X_{2}$ (or $X_{2}$ given $X_{1}$ respectively). In the following we will do both.

For the marginal density of $X_{1}$ we need to obtain $f_{X_1}(x_{1}) = \int_{\mathbb{R}} f(x_{1}, x_{2}) \mbox{d}x_{2}$ and, analogously, for the marginal density of $X_{2}$, one needs to integrate the join density wrt. the first covariate instead, i.d., $f_{X_2}(x_{2}) = \int_{\mathbb{R}} f(x_{1}, x_{2}) \mbox{d}x_{1}$. Both marginals are again normaly distributed and it holds that $X_{1} \sim N(\mu_1, \sigma_1^2)~~~$ and $~~~X_{2} \sim N(\mu_2, \sigma_2^2)$.
For a simple example with a two dimensional normal distribution the conditional distribution distribution of $X_{2}$ given $X_{1} = x_{1}$ is, again, normal and it holds that \[(X_{2} | X_{1 } = x_{1}) \sim N\Big(\mu_{2} + \frac{\sigma_{21}(x_1 - \mu_1)}{\sigma_{1}^2}, \sigma_{2}^2 - \frac{\sigma_{12}\sigma_{21}}{\sigma_{1}^2}\Big).\] An analogous expression can be obtained also for the conditional distribution of $X_{1}$ given $X_{2}$.

Now we can apply the formulas given above to obtain the marginal and conditional distributions for some given two-dimensional normal distribution. We will use the R library mvtnorm (which needs to be firstly installed on R). The library is loaded into the working environment by running the command

library("mvtnorm")

Let us consider a simple example with two dimensional normal distribution with the zero mean vector $\boldsymbol{\mu} = (0,0)^\top$, and the variance-covariance matrix $\Sigma = \left( \begin{array}{cc} 1 & 0.8 \\0.8& 1\end{array} \right)$. We would like to calculate the conditional distribution of $X_{2}$ given $X_{1} = 0.7$.

Individual task

Is there any linear relationship between the covariates $X_1$ and $X_2$? Can you quantitatively express how strong this relationship is?

In terms of the linear regression modelling approach: imagine you obtain a sample from the given two dimensional normal distribution and you fit a simple regression line to the data. Do you have some expectation about the parameter estimates you obtain when fitting the linear regression model?
What about the correlation coefficient (respectively its sample version) and the Multiple $R^2$ value in the summary output?

n <- 100
sample <- rmvnorm(n, c(0, 0), matrix(c(1, 0.8, 0.8, 1),2,2))
(summary(lm(sample[,1] ~ sample[,2])))

## 
## Call:
## lm(formula = sample[, 1] ~ sample[, 2])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.98635 -0.33903 -0.03215  0.29492  1.51192 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.01122    0.04849   0.231    0.818    
## sample[, 2]  0.81707    0.04474  18.263   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4789 on 98 degrees of freedom
## Multiple R-squared:  0.7729, Adjusted R-squared:  0.7706 
## F-statistic: 333.5 on 1 and 98 DF,  p-value: < 2.2e-16

(cor(sample[,1], sample[,2]) * sd(sample[,1]) /sd(sample[,2]))

## [1] 0.8170716

plot(sample[,1] ~ sample[,2], pch = 21, bg = "gray")
abline(lm(sample[,1] ~ sample[,2]), col = "red", lwd = 2)
lines(c(-4, 4) ~ rep(mean(sample[,2]), 2), col = "blue", lty = 2)
lines(rep(mean(sample[,1]), 2) ~ c(-4, 4), col = "blue", lty = 2)

Note, that the blue dashed lines do not represent the point $[0,0] \in \mathbb{R}^2$. Instead, they denote c(mean(sample[,2]), mean(sample[,1])).

Use the following piece of the R code to obtain a comparison between the joint distribution, marginal distribution, and the conditional distribution of $X_{2}$ given $X_{1} = 0.7$. Derive the theoretical expressions for both, the marginal distributions and the conditional distribution.

Sigma <- matrix(c(1,.8,.8,1), nrow=2) ## variance-covariance matrix

x <- seq(-3,3,0.01)
contour(x,x,outer(x,x,function(x,y){dmvnorm(cbind(x,y),sigma=Sigma)}), col = "blue")

abline(v=.7, lwd=2, lty=2, col = "red")
text(0.75, -2, labels=expression(x[1]==0.7), col = "red", pos = 4)

### conditional distribution of X2 | X1 = 0.7
y <- dnorm(x, mean =  0.8 * 0.7, sd = sqrt(1 - 0.8^2))
lines(y-abs(min(x)),x,lty=2,lwd=2, col = "red")

### marginals
m1 <- m2 <- dnorm(x, 0, 1)
lines(x, m1 - abs(min(x)), lty = 1, lwd = 2, col = "gray30")
lines(m2 - abs(min(x)), x, lty = 1, lwd = 2, col = "gray30")

The conditional distribution can be obtained for any value $X_{1} = x_1$, for instance, in the following we vizualize the conditional distribution of $(X_{2} | X_{1} = -1)$. Marginal distributions remain both the same. The conditional distribution changed correspondingly.

contour(x,x,outer(x,x,function(x,y){dmvnorm(cbind(x,y),sigma=Sigma)}), col = "blue")
abline(v=-1, lwd=2, lty=2, col = "red")

### marginals
lines(x, m1 - abs(min(x)), lty = 1, lwd = 2, col = "gray30")
lines(m2 - abs(min(x)), x, lty = 1, lwd = 2, col = "gray30")

### conditional distribution given X1 = 0.7
abline(v=.7, lwd=2, lty=2, col = "gray")
lines(y + 0.7,x,lty=2,lwd=2, col = "gray")
lines(y - abs(min(x)),x,lty=2,lwd=2, col = "gray")
text(0.75, -2, labels=expression(x[1]==0.7), col = "gray", pos = 4)

### conditional distribution of X2 | X1 = - 1
y2 <- dnorm(x, mean = 0.8 * (- 1), sd = sqrt(1 - 0.8^2))
lines(y2 - abs(min(x)),x,lty=2,lwd=2, col = "red")

Another option in R to obtain the same output:

contour(x,x,outer(x,x,function(x,y){dmvnorm(cbind(x,y),sigma=Sigma)}), col = "blue")
abline(v=-1, lwd=1, lty=1, col = "red")
lines(y2 - 1 ,x,lty=2,lwd=2, col = "red")

And even more generalized version for a set of different values of $x_{1} \in \mathbb{R}$:

contour(x,x,outer(x,x,function(x,y){dmvnorm(cbind(x,y),sigma=Sigma)}), col= "blue")
condN <- function(x, cx){dnorm(x, mean = 0.8 * cx, sd = sqrt(1 - 0.8^2))}
for (i in seq(-2, 2, by = 0.25)){col <- colors()[grep("red",colors())][4*i + 9]; 
                                 lines(condN(x, i) + i, x, lwd = 2, col = col);
                                 abline(v = i, lwd=1, lty=1, col = col)}

The theoretical model above can be now compared with the sample version below:

set.seed(1234)
s <- rmvnorm(5000, c(0, 0), Sigma)
plot(s, pch = 21, bg = "lightblue", xlim = c(-3, 3), ylim = c(-3, 3), xlab = "", ylab = "", cex = 0.6)
abline(lm(s[,2] ~ s[,1]), col = "red", lwd = 2)

set.seed(1234)
s <- rmvnorm(5000, c(0, 0), Sigma)
plot(s, pch = 21, bg = "lightblue", xlim = c(-3, 3), ylim = c(-3, 3), xlab = "", ylab = "", cex = 0.8)
abline(lm(s[,2] ~ s[,1]), col = "red", lwd = 2)

for (i in seq(-2, 2, by = 0.25)){col <- colors()[grep("red",colors())][4*i + 9]; 
                                 lines(condN(x, i) + i, x, lwd = 2, col = col);
                                 abline(v = i, lwd=1, lty=1, col = col)}

Individual task

Consider the conditional distribution we obtained above. Can you numerically reconstruct (obtain) the values which are used as the parameters in the normal distribution?

Recall the general formula for a conditional distribution of a random vector $\boldsymbol{X}_{2}$ given $\boldsymbol{X}_{1} = \boldsymbol{x}_{1}$ if the join distribution of $(\boldsymbol{X}_{1}^\top, \boldsymbol{X}_{2}^\top)^\top$ is normal and $\boldsymbol{X}_{1}$ is a $p$-dimensional random vector for $p > 1$ and $\boldsymbol{X}_{2}$ is a $q$-dimensional random vector with $q > 1$.

Recal a normal linear model $Y_{i} = \alpha + \beta X_i + \varepsilon_i$, for $i = 1, \dots, n$. What is the distribution of $Y$ given $X = x$? Do you see the relationship of the model with the countour plot above?

Consider a multivariate example for at least a four dimensional normal distribution and calculate the conditional distribution of the first two elements $(X_{1}, X_{2})$ given the remaining elements. Create a contour plot for the conditional two-dimensional normal distribution with various values for the condition.
Recall, that the joint distribution of some random vector (for instance, a normally distributed random vector) can be straightforwardly used to obtain any marginal distribution of interest. Analogously, the corresponding conditional distributions can be formed. On the other hand, however, having all possible marginal distributions is still not enough to reconstruct the whole joint distribution. Specifically, normally distributed marginals of some random vector do not imply normality of the whole random vector.

Comment

For a general multivariate normal distribution denoted as

$\Big(\begin{array}{x}\boldsymbol{X}_{1}\\\boldsymbol{X}_{2}\end{array}\Big) \sim N_{p + q}\left(\boldsymbol{\mu} = \Big(\begin{array}{c} \boldsymbol{\mu}_{1} \\ \boldsymbol{\mu}_{2}\end{array}\Big), \Sigma = \left( \begin{array}{cc} \Sigma_{11} & \Sigma_{21} \\\Sigma_{12} & \Sigma_{22} \end{array} \right) \right)$

for some random vectors $\boldsymbol{X}_{1} \in \mathbb{R}^{p}$ and $\boldsymbol{X}_{2} \in \mathbb{R}^{q}$ for $p, q \in \mathbb{N}$ and the corresponding expectation vectors $\boldsymbol{\mu}_{1}$ and $\boldsymbol{\mu}_{2}$ and the variance-covariance matrix $\Sigma$ the conditional distribution of $\boldsymbol{X}_{1}$ given $\boldsymbol{X}_{2} = \boldsymbol{x}_{2}$ is defined by the following expression:

$(\boldsymbol{X}_{1} | \boldsymbol{X}_{2} = \boldsymbol{x}_{2}) \sim N_{p}\left(\boldsymbol{\mu}_{1} + \Sigma_{12}\Sigma_{22}^{-1}(\boldsymbol{x}_{2} - \boldsymbol{\mu}_{2} ), \Sigma = \Sigma_{11} - \Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21} \right)$.

2. Conditional multivariate normal distribution with the ‘condMVNorm()’ package

There is a nice (and quite brief) R library available from the R cran repository which can be used to calculate the conditional distribution of some random vector $\boldsymbol{X}_{1}$ given another random vector (value) $\boldsymbol{X}_{2} = \boldsymbol{x}_{2}$, where the joint distribution of $(\boldsymbol{X}_{1}^{\top}, \boldsymbol{X}_{2}^{\top})^{\top}$ is assuned ti be multivariate normal (with the corresponding dimension and the given parameters). The library can be installed by using command install.packages("condMVNorm").

library("condMVNorm")

Functions dcmvnorm() and condMVN can be directly used to obtain the conditional distribution of interest (see the help session and use ?dcmvnorm and ?condMVN for more details).

We apply these functions to the previous example with a two dimensional normal distribution with a zero mean vector and the variance-covariance matrix $\Sigma = \left( \begin{array}{cc} 1 & 0.8 \\ 0.8 & 1 \end{array} \right)$. We want to obtain the conditional distribution of $X_{2}$ given $X_1 = 0.7$.

condDist <- condMVN(mean = c(0,0), sigma = Sigma, dependent.ind = 2, given.ind = 1, X.given = 0.7)

and also

dcmvnorm(0, mean = c(0,0), sigma = Sigma, dependent.ind = 2, given.ind = 1, X.given = 0.7)

## [1] 0.4301297

which should be compared with

dnorm(0, condDist$condMean, sqrt(condDist$condVar))

## [1] 0.4301297

Comment

Keep in minds that some functions in R require standard deviation to be specified when working with the normal distribution (e.g. functions dnorm(), pnorm(), qnorm(), etc.) while others need the variance-covariance matrix (variance in one dimension) to be provided instead.

3. Conditional normal distribution – illustration example

In this example we again consider the well-known dataset mtcars which is standardly available in the R environment (for a more detailed description of the dataset see the help session ?mtcars).

summary(mtcars)

##       mpg             cyl             disp             hp       
##  Min.   :10.40   Min.   :4.000   Min.   : 71.1   Min.   : 52.0  
##  1st Qu.:15.43   1st Qu.:4.000   1st Qu.:120.8   1st Qu.: 96.5  
##  Median :19.20   Median :6.000   Median :196.3   Median :123.0  
##  Mean   :20.09   Mean   :6.188   Mean   :230.7   Mean   :146.7  
##  3rd Qu.:22.80   3rd Qu.:8.000   3rd Qu.:326.0   3rd Qu.:180.0  
##  Max.   :33.90   Max.   :8.000   Max.   :472.0   Max.   :335.0  
##       drat             wt             qsec             vs        
##  Min.   :2.760   Min.   :1.513   Min.   :14.50   Min.   :0.0000  
##  1st Qu.:3.080   1st Qu.:2.581   1st Qu.:16.89   1st Qu.:0.0000  
##  Median :3.695   Median :3.325   Median :17.71   Median :0.0000  
##  Mean   :3.597   Mean   :3.217   Mean   :17.85   Mean   :0.4375  
##  3rd Qu.:3.920   3rd Qu.:3.610   3rd Qu.:18.90   3rd Qu.:1.0000  
##  Max.   :4.930   Max.   :5.424   Max.   :22.90   Max.   :1.0000  
##        am              gear            carb      
##  Min.   :0.0000   Min.   :3.000   Min.   :1.000  
##  1st Qu.:0.0000   1st Qu.:3.000   1st Qu.:2.000  
##  Median :0.0000   Median :4.000   Median :2.000  
##  Mean   :0.4062   Mean   :3.688   Mean   :2.812  
##  3rd Qu.:1.0000   3rd Qu.:4.000   3rd Qu.:4.000  
##  Max.   :1.0000   Max.   :5.000   Max.   :8.000

attach(mtcars)

For illustrative purposes only, we will simply assume that (continuous) covariates ‘mpg’ (miles per gallon), ‘hp’ (horse power) and ‘disp’ (displacement) are all normaly distributed—the triple $(X_{1}, X_{2}, X_{3})^{\top} = (mpg, hp, disp)^{\top}$ follows a three dimensional normal distribution with some unknown mean vector $\boldsymbol{\mu} = (\mu_{1}, \mu_{2}, \mu_{3})^{\top} \in \mathbb{R}^{3}$ and some positive definite and symmetric variance-covariance matrix $\Sigma$.

Use the R software to obtain the corresponding estimates for the mean vector $\boldsymbol{\mu} = (\mu_{1}, \mu_{2}, \mu_{3})^{\top} \in \mathbb{R}^{3}$ and the variance-covariance matrix $\Sigma$.

What is the (estimated) distribution if the miles per gallon (‘mpg’) given the horse power and displacement, where $(X_{2}, X_{3})^{\top} = (120,120)^{\top}$?

The corresponding estimates can be obtained by

mu <- apply(cbind(mpg, hp, disp), 2, mean) 
Sigma <- cov(cbind(mpg, hp, disp))

Thus, the corresponding conditional mean and variance are obtained as

condMean <- mu[1] + Sigma[1, 2:3] %*% solve(Sigma[2:3,2:3]) %*% (c(120, 120) - mu[2:3])
condSigma <- Sigma[1,1] - Sigma[1, 2:3] %*% solve(Sigma[2:3,2:3]) %*% Sigma[2:3, 1]

and given the assumptions of multivariate normality for the triple $(X_{1}, X_{2}, X_{3})^{\top} = (mpg, hp, disp)^{\top}$ we get the univariate normal distribution with the following moments:

condMean

##          [,1]
## [1,] 24.11354

condSigma

##         [,1]
## [1,] 9.14495

plot(dnorm(x <- seq(0, 35, length = 500), mean = condMean, sd = sqrt(condSigma)) ~ x, 
                                           xlab = "Miles per gallon", ylab = "Conditional Density", 
                                           type = "l", lwd = 2, col = "red", 
                                           xlim = c(0, 35), ylim = c(0,0.18))
lines(dnorm(x, mean(mpg), sd(mpg)) ~ x, lwd = 1)
legend(0, 0.175, legend = c("Conditional distribution: MPG | HP = 120 & Disp = 120", 
                           "Estimated Normal Distribution of MPG"), col = c("red", "black"), lty = c(1,1))

Individual task

Consider a simple linear regression model where the miles per gallon information is used as a dependent covariate and the horse power and the displacement value are both used as the explanatory variables. Estimate the model and interpret the estimated coefficient.
```
m <- lm(mpg ~ hp + disp, data = mtcars)
```
Use the model above and estimate the expected value of the consuption (miles per gallon) for a car with the horse power being equal to 120 and the displacement being 120 too.
```
CondMean <- m$coeff[1] + m$coeff[2] * 120 + m$coeff[3] * 120
```
Compare the conditional distribution you obtained above with the prediction obtained by the linear regression model. Discuss the results.

Theoretical Exercises

Exercise 1
Let $\boldsymbol{X}_1 \sim N_{q}(\boldsymbol{\mu}_1, \boldsymbol{\Sigma}_{11})$ and $(\boldsymbol{X}_2 | \boldsymbol{X}_1 = \boldsymbol{x}_1) \sim N_{p - q}(\mathbb{A}\boldsymbol{x}_1 + \boldsymbol{b}, \boldsymbol{\Omega})$, where $\boldsymbol{\Omega}$ does not depend on $\boldsymbol{x}_1 \in \mathbb{R}^q$. Show, that the joint distribution of the whole $p$-dimensional random vector $(\boldsymbol{X}_1^\top, \boldsymbol{X}_2^\top)^\top$ satisfy \[ (\boldsymbol{X}_1^\top, \boldsymbol{X}_2^\top)^\top \sim N_{p}(\boldsymbol{\mu}, \boldsymbol{\Sigma}), \] where \[ \boldsymbol{\mu} = \left( \begin{array}{c} \boldsymbol{\mu}_1\\ \mathbb{A}\boldsymbol{\mu}_1 + \boldsymbol{b} \end{array} \right) \qquad \textrm{and} \qquad \boldsymbol{\Sigma} = \left( \begin{array}{cc} \boldsymbol{\Sigma}_{11} & \boldsymbol{\Sigma}_{11} \mathbb{A}^\top\\ \mathbb{A}\boldsymbol{\Sigma}_{11} & \boldsymbol{\Omega} + \mathbb{A}\boldsymbol{\Sigma}_{11}\mathbb{A}^\top \end{array} \right). \]

Exercise 2
Consider the random vector $\boldsymbol{X}$ with the two dimensional normal distribution $N_{2} \left( \left( \begin{array}{c}1\\2\end{array} \right), \left( \begin{array}{cc} 2 & 1\\1 & 2\end{array} \right) \right)$ and the conditional random vector $\boldsymbol{Y} | \boldsymbol{X}$, such that \[ \boldsymbol{Y} | \boldsymbol{X} \sim N_{2}\left( \left( \begin{array}{c} X_{1}\\ X_{1} + X_{2} \end{array} \right), \left( \begin{array}{cc} 1 & 0\\ 0 & 1 \end{array} \right) \right). \]
Determine the distribution of $Y_2 | Y_{1}$ and the distribution of $\boldsymbol{W} = \boldsymbol{X} - \boldsymbol{Y}$.

Exercise 3
Consider a random vector $\boldsymbol{X} \sim N_{2}(\boldsymbol{\mu}, \Sigma)$, for $\boldsymbol{\mu} = (2,2)^\top$ and $\Sigma = \mathbb{I}_2$ (a unit matrix). Consider also matrices $\mathbb{A} = (1, 1)$ and $\mathbb{B} = (1, -1)$. Show, that the random variables $\mathbb{A}\boldsymbol{X}$ and $\mathbb{B}\boldsymbol{X}$ are independent.

Exercise 4
Consider a random vector $(X, Y, Z)^\top \sim N_3(\boldsymbol{\mu}, \Sigma)$. Compute the values for $\boldsymbol{\mu}$ and $\Sigma$, if you know, that
- $Y | Z \sim N(-Z, 1)$;
- $\mu_{Z|Y} = -\frac{1}{3} - \frac{1}{3}Y$;
- $X | (Y, Z) \sim N(2 + 2Y + 3Z, 1)$;
Determine the conditional distribution of $X|Y$ and also $X | Y + Z$.

Exercise 5
Assume, you know the following:
- $Z \sim N(0,1)$;
- $Y | Z \sim N(1 + Z, 1)$;
- $X | (Y, Z) \sim N(1 - Y, 1)$;
Find the distribution of the random vector $(X, Y, Z)^\top$ and the conditional distribution of $Y | (X, Z)$. Consider the linear transformations $U = 1 + Z$ and $V = 1 - Y$ and find the joint distribution of the random vector $(U, V)^\top$. Calculate the value of $E(Y | U = 2)$.

Homework Assignment

(Deadline: Lab session no. 4 / 08.03.2022)

Choose one theoretical example from the list above and solve it. Use the R simulations to empirically verify the results. Graphical outputs of the simulations are to be provided on the web site.
As already mentioned above, the joint normality of some random vector directly implies normality of all possible marginals. On the other hand, the normal distribution of all marginals does not necessarily imply the joint normality of the whole random vector.

Try to propose a simple example (using a two/three dimensional random vector $(X_1, X_2)^\top$ or $(X_1, X_2, X_3)^\top$ respectively) showing that despite all normaly distributed marginals the whole random vector fails to follow the joint normal distribution.
Think of some appropriate options to vizualize the whole problem. Graphical outputs are again to be provided on the web site.
Provide the solution on you personal web page by Thuesday, 08.03.2022, 12:00 at latest.