<table width = 900 border = 0><tr><td>
# <b>NMST539 | Lab Session 3</b>

## <b>Multivariate Normal Distribution</b>
### (marginal and conditional distributions)
### <b><font color = "#880000">LS 2017 | Monday 05/03/18</font></b>
###### <a href = "cvicenie2018_3.Rmd">Rmd file (UTF8 coding)</a>

**************

The R-software is available for download from the website: <a href = "https://www.r-project.org">https://www.r-project.org</a>

A user-friendly interface (one of many): <a href = "https://www.rstudio.com">RStudio</a>.<br><br>

<b>Manuals and introduction into <i>R</i>
(in Czech or English):</b><br>
		
<ul>
<li>Bína, V., Komárek, A. a Komárková, L.: <i>Jak na jazyk R.</i>  (<a href = "NMFM301/Rmanual2.pdf">PDF súbor</a>)</li>
<li>Komárek, A.: <i>Základy práce s R.</i>  (<a href = "http://msekce.karlin.mff.cuni.cz/~komarek/vyuka/dataRko/Rmanual1.pdf">PDF súbor</a>)</li>
<li>Kulich, M.: <i>Velmi stručný úvod do R.</i> (<a href = "NMFM301/uvodr.pdf" >PDF súbor</a>)</li>
<li>De Vries, A. a Meys, J.: <i>R for Dummies.</i> (ISBN-13: 978-1119055808) </li>		
</ul>
**************

#### <font color = "#880000"><b>1. Conditional Normal Distribution</b></font>

Let us consider a two-dimensional normal distribution of some random vector $\Big(\begin{array}{x}X_{1}\\X_{2}\end{array}\Big)$. The corresponding distribution is usually  denoted as 

<center>
$\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)$,
</center>

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

<center>
$\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\},}$
</center><br>

for an arbitrary $\boldsymbol{x} = (x_{1}, x_{2})^{\top} \in \mathbb{R}^{2}$.<br><br>

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. 

<ul>
<li>For the marginal density of $X_{1}$ we need to obtain  $f(x_{1}) = \int_{\mathbb{R}} f(x_{1}, x_{2}) \mbox{d}x_{2}$ and analogously also for the marginal density of $X_{2}$, where  integrate the join density wrt the first covariate instead. Both marginals are again normaly distributed and it holds that 
<center>
$X_{1} \sim N(\mu_1, \sigma_1^2)~~~$ and $~~~X_{2} \sim N(\mu_2, \sigma_2^2)$.
</center><br><br>
</li>

<li>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 (analogously also for the distribution of $X_{1}$ given $X_{2}$)

<center>
$(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).$
</center></li>
</ul>


<br><br>




Now we can apply the formulas given above to obtain the marginal and conditional distributions. 
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 

```{r, message=F, warning=F}
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$.<br><br>

#### <font color = "#880000">Do by Yourselves</font>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<ul>
<li>Is there any linear relationship between the covariates $X_1$ and $X_2$? Can you quantitatively express how strong this relationship is?</li><br>

<li>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? <br>
Try the following:

```{r, eval = F}
n <- 100
sample <- rmvnorm(n, c(0, 0), matrix(c(1, 0.8, 0.8, 1),2,2))
summary(lm(sample[,1] ~ sample[,2]))
```
</li>
</ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<br><br>


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 the marginal distributions and the conditional distribution.  


```{r, fig.width = 10}
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, we obtain the conditional distribution of $X_{2} | X_{1} = -1$): 

```{r, fig.width = 10}
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")

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

```

Another option in R to obtain the same output:  

```{r, fig.width = 10}
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}$:

```{r, fig.width = 10}
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)}
```

<br><br>

#### <font color = "#880000">To Do by Yourself</font>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<ul>
<li>Consider the conditional distribution we obtained above. Can you numerically reconstruct (obtain) the values which are used as the parameters in the normal distribution?</li><br>

<li>Recall the general formula for a conditional distribution of $X_{2}$ given $X_{1} = x_{2}$ if the join distribution of $(X_{1}^\top, X_{2}^\top)^\top$ is normal 
and $X_{1}$ is $p$-dimensional random vector for $p > 1$ and $X_{2}$ is $q$-dimensional random vector with $q > 1$.</li><br>

<li>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? </li><br>

<li>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.</li>
</ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />

<br>

#### <font color = "#880000">Comment</font>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
For a general multivariate normal distribution denoted as 

<center>
$\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)$
</center><br><br>

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: 

<center>
$(\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)$.
</center>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<br><br>

#### <font color = "#880000"><b>2. Conditional Multivariate Normal with 'condMVNorm()' package</b></font>

There is a nice (and quite brief) R library available from the R cran mirrors 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 multivariate normal (with the corresponding dimension). The library can be installed by using command `install.packages("condMVNorm")`. <br><br>

```{r}
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$.

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

and also

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

which should be compared with 

```{r}
dnorm(0, condDist$condMean, sqrt(condDist$condVar))
```

<br>

#### <font color = "#880000">Comment</font>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
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. 

<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<br>


#### <font color = "#880000"><b>3. Conditional Normal Distribution Example</b></font>

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

```{r}
summary(mtcars)
attach(mtcars)
```

Let us assume that covariates 'mpg' (miles per gallon), 'hp' (horse power) and 'disp' (displacement)
are all normaly distributed and 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$. 

<ul>
<li>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$.</li>
</ul>

```{r, echo=FALSE, fig.width = 10}
par(mfrow = c(1,3))
hist(mpg, xlab = "Miles per gallon", ylab = "Frequency", col = "lightblue")
hist(hp, xlab = "Gross horsepower", ylab = "", col = "lightblue")
hist(disp, xlab = expression(paste("Displacement [", cu^3, "]", sep = "")), ylab = "", col = "lightblue")
```

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 

```{r}
mu <- apply(cbind(mpg, hp, disp), 2, mean) 
Sigma <- cov(cbind(mpg, hp, disp))
```

Thus, the corresponding conditional mean and variance are obtained as

```{r}
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]
```

Which gives us the following distribution function: 

```{r, fig.width = 10}
condMean
condSigma

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))
```


#### <font color = "#880000">To Do by Yourself</font>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<ul>
<li>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 explanatory variables. Estimate the model and interpret the estimated coefficient.</li>

<li>Estimate the expected value of the consuption (miles per gallon) for a car with the horse power 120 and the displacement 120 too. </li>

<li>Compare the conditional distribution you obtained above with the prediction obtained by the linear regression model. Discuss the results.</li>

</ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />

<br><br>

#### <font color = "#880000">Theoretical Exercises</font>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<ul>
<li><b>Exercise 1</b><br> 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}$.</li><br><br>

<hr size="3" width="100%" noshade style="color:#440000" align="left" />

<li><b>Exercise 2</b><br> 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)^\top$ and $\mathbb{B} = (1, -1)^\top$. Show, that the random variables $\mathbb{A}\boldsymbol{X}$ and $\mathbb{B}\boldsymbol{X}$ are independent.</li><br><br>

<hr size="3" width="100%" noshade style="color:#440000" align="left" />

<li><b>Exercise 3</b><br> 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 

<ul>
<li>$Y | Z \sim N(-Z, 1)$;</li>
<li>$\mu_{Z|Y} = -\frac{1}{3} - \frac{1}{3}Y$;</li>
<li>$X | (Y, Z) \sim N(2 + 2Y + 3Z, 1)$;</li>
</ul>
Determine the conditional distribution of $X|Y$ and also $X | Y + Z$.
</li><br><br>

<hr size="3" width="100%" noshade style="color:#440000" align="left" />

<li><b>Exercise 4</b><br> Assume, you know the following: 
<ul>
<li>$Z \sim N(0,1)$;</li>
<li>$Y | Z \sim N(1 + Z, 1)$;</li>
<li>$X | (Y, Z) \sim N(1 - Y, 1)$;</li>
</ul>
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)$.
</li>

</ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />




<br><br>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<hr size="3" width="100%" noshade style="color:#440000" align="left" /><br><br> 


## <font color = "#880000"><b>Homework Assignment</b></font>
### (Deadline: lab session no. 4 / 12.03.2018)

Use the command below to instal the `SMSdata` packages with various multivariate datasets: 

```{r, eval = F}
install.packages(pkgs="http://www.karlin.mff.cuni.cz/~hlavka/sms2/SMSdata_1.0.tar.gz", repos=NULL, type="source")
```

Chose one dataset of your choise from the list of all available datasets in the package:

```{r, eval = F}
data(package = "SMSdata") 
```
 
 There are 21 different datasets and you can load each of them by typing and running 
`data(dataset_name)`, for instance:

```{r, eval = F}
data(plasma) 
```

<ul>
<li>Choose one dataset where there are at least two covariates which can be considered to be continous ones wich (approximatelly) some normal distribution (e.g. use a histogram 
plot to help yourself to decide) and load the data into the R environment. </li>

<li>Use your `Rnw` file (which you have already prepared for the previous homework assignment) and add some figures into your report. At least one figure should contours of the two-dimensional distribution you are considering. Assuming the two-variate join normal distribution for your two covariates calculate also a conditional distribution when conditioning one covariate by another one and both marginal distribtions. Draw the corresponding distributions into the firgure.</li>

<li>Provide some reasonable interpretation of the figures and all elements included in them.</li>
</ul>



</td></tr></table>
<br><br><br>