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

## <b>Multivariate Normal Distribution</b>
### <font color = "#880000">LS 2017 | Thursday 02/03/17</font>
###### <a href = "cvicenie2017_2.Rmd">Rmd file (UTF8 coding)</a>

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

For the second NMST539 lab session we will discover some tools, available within the R software instalation, which are meant
for working with the multivariate normal distribution. 

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><br>

#### <font color = "#880000"><b>R package 'mvtnorm'</b></font>

For the beginning, we need to make sure that the library for the multivariate normal distribution (`mvtnorm`) is loaded in the R working environment (use commad `library("mvtnorm")`). If the library can not be loaded we need to install it first by running the command `install.packages("mvtnorm")`

```{r}
library("mvtnorm")
```

We can now use the loaded library to easily generate random values from some multivariate normal distribution with some pre-specified variance-covariance matrix $\Sigma$. 

For more details about the library (including the manual document) see the library website on R cran:<br> <a href = "https://cran.r-project.org/web/packages/mvtnorm/index.html">https://cran.r-project.org/web/packages/mvtnorm/index.html</a>
   
<br><br>   
   
#### <font color = "#880000">Comment</font>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
It is useful to always keep in mind that everytime we are about to use some random generator in R (or any other software) it is good to predefine the initial values for the generator (in R there is a command called `set.seed()` to do the job). By doing so one makes sure that the obtained  results are always all reconstructable and they can be easily verified. 

For instance try (among each other) the command:

```{r}
(sample <- rnorm(10))
```

and compare it with the same commend, however, with the initial setting of the generator done by the `set.seed(1234)` command (any numerical value can be used as a parameter): 

```{r}
set.seed(1234)
(sample <- rnorm(10))
```

Is it clear, what is the difference between these two outputs? 
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<hr size="3" width="100%" noshade style="color:#440000" align="left" /><br><br>   

For the beginning, let us start with a simple (two-dimensional only) example: we will generate a random sample of size $n \in \mathbb{N}$ from the two-dimensional normal distribution

$N_{2}\Big(\boldsymbol{\mu} = \Big(\begin{array}{c}2\\3\end{array}\Big), \Sigma = \Big(\begin{array}{cc}10^2 & 6^2\\ 6^2 & 6^2\end{array}\Big) \Big)$.

Remember, that in case of a one-dimensional random generator for the normal distribution in R (command `rnorm()`) one needs to specify the standard error instead of the variance (compare the help sessions `?rnorm` and `?rmnorm`). 

The sample of size $n=1000$ from the given two-dimensional normal distribution can be generated by the following command: 

```{r}
set.seed(1234)
s1 <- rmvnorm(1000, c(2, 3), matrix(c(10^2, 6^2, 6^2, 6^2),2,2))
```
 
And we can plot the generated random sample into a nice scatterplot using for instance the command

```{r, fig.width = 10}
plot(s1, pch=21, xlim=c(-35, 35), ylim=c(-35,35), xlab="Marginal No. 1", ylab="Marginal No. 2", bg = "lightblue")
```

Using now mean estimates and the sample variance-covariance matrix we can easily draw some ‘depth’ contours in the figure. Another option of course, is to use the theoretical values, as we know the distribution we generated from (which is usually not the case in real life situations). However, in the following, we will rather use the sample  and we calculate the sample mean vector and the sample variance-covariance matrix.

```{r}
center <- apply(s1, 2, mean)   ## sample mean vector
sigma <- cov(s1)               ## sample variance-covariance matrix
```

The two dimensional normal density function is given by

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

where $\boldsymbol{x} = (x_{1}, x_{2})^\top \in \mathbb{R}^2$ and $\boldsymbol{\mu} = (\mu_{1}, \mu_{2})^\top$ is the mean vector.

Thus, we also need the inverse matrix for the sample variance-covariance matrix $\widehat{\Sigma}$ to be able to draw the corresponding contours. The R fuction `ellipse()` in the following code is responsible for calculating the value of the two dimensional normal density function when plugging in the appropriate sample estimates and the margin values. 

```{r}
sigma.inv = solve(sigma, matrix(c(1,0,0,1),2,2))

ellipse <- function(s,t){
                    u<-c(s,t)-center; 
                    e <- (-1) * (u %*% sigma.inv %*% u) / 2; 
                    exp(e)/(2 * pi * sqrt(det(sigma.inv)))}
```

This can be used as follows: 

```{r, fig.width = 10}
n <- 60
x <- (0:(n-1))*2 - 50
y <- (0:(n-1))*2 - 50

z <- mapply(ellipse, as.vector(rep(x,n)), as.vector(outer(rep(0,n), y, `+`)))
```

And the final plot is produced by the following couple of the commands:

```{r, fig.width = 10}
plot(s1, pch=21, xlim=c(-35, 35), ylim=c(-35,35), xlab="Margin No. 1", ylab="Margin No. 2", bg = "lightgray")
contour(x,y,matrix(z,n,n), levels=(0:15), col = terrain.colors(16), add=TRUE, lwd = 1)
```


<br><br>

#### <font color = "#880000">Do by Yourselves</font>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<ul>
<li>Try various settings for the parameters defining the two dimensional normal distribution (the mean vector and the variance covariance matrix) used above and try to understand what happens with the distribution of the points, rotation of the structure and the corresponding countours when changing the corresponding parameters.</li><br>

<li>Use two different samples each simulated from the two dimensional normal distribution with different parameters and draw the samples with the corresponding countours into one plot only. Again try to use different parameter settings to see what happens.</li><br>

<li>Try to give a resonable and sensible interpretation (in a way that non-statisticians or non-mathemticians would understant) of the parameters used in the models.</li>
</ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />

<br><br>

Another option which can be used to get the same results is the function predefined in the R package `mixtools`. Install the library by using `install.packages(mixtools)` and load the library by running the command 

```{r, warning=F, message=F}
library("mixtools")
```

Within the `mixtools` library there is another `ellipse()` function defined (use the help session `?ellipse` to see how it works). 

```{r, fig.width = 10}
rm(list = ls())
library(mixtools)
library(mvtnorm) 

set.seed(1234)
p <- rmvnorm(1000, c(250000, 20000), matrix(c(100000^2, 22000^2, 22000^2, 6000^2),2,2))
plot(p, pch=20, xlim=c(0,500000), ylim=c(0,50000), xlab="Margin No. 1", ylab="Margin No. 2")
ellipse(mu=colMeans(p), sigma=cov(p), alpha = .05, npoints = 250, col="red") 
```

<br><br>

#### <font color = "#880000">Do by Yourselves</font>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
Can you reconstruct the countour from the figure above? How is the red ellipse calculated there?

```{r, fig.width = 10}
alpha <- 0.05
angles <- seq(0, 2*pi, length = 200)
cpoints <- cbind(cos(angles), sin(angles))

q <- sqrt(qchisq(1 - alpha, df = 2))
mu <- colMeans(p)
sigma <- cov(p)
sigmaEig <- eigen(sigma)
sigmaSqr <- sigmaEig$vectors %*% diag(sqrt(sigmaEig$values)) %*% t(sigmaEig$vectors)
Epoints <-  q * (cpoints %*% sigmaSqr)
Epoints[,1] <- Epoints[,1] + mu[1]
Epoints[,2] <- Epoints[,2] + mu[2]

plot(p, pch=20, xlim=c(0,500000), ylim=c(0,50000), xlab="Margin No.1", ylab="Margin No.2")
points(Epoints[,1], Epoints[,2], col = "red", pch = "*")
```
<hr size="3" width="100%" noshade style="color:#440000" align="left" />

<br><br>

#### <font color = "#880000"><b>Alternative Approaches with the library `mvrnorm()`</b></font>

Another option for simulating a general multivariate normal sample is available within the library `MASS` by running the command `mvrnorm()` (again, use the help session `?mvrnorm` to get into more details). The library needs to be loaded in by the standard command

```{r, message=F, warning=FALSE}
library("MASS")
```

Next, we can generate some random samples from two-dimensional normal distribution and we can even calculate the corresponding 
multivariate kernel density estimates. 

```{r, fig.width=10}
### two independent samples from two-dimensional normal each
bivn1 <- mvrnorm(1000, mu = c(0, 0), Sigma = matrix(c(1, .5, .5, 1), 2))
bivn2 <- mvrnorm(1000, mu = c(0, 0), Sigma = matrix(c(1.5, 1.5, 1.5, 1.5), 2))

### two dimensional kernel density estimate
bivn.kde1 <- kde2d(bivn1[,1], bivn1[,2], n = 50)
bivn.kde2 <- kde2d(bivn2[,1], bivn2[,2], n = 50)

par(mfrow = c(1,2))
persp(bivn.kde1, phi = 45, theta = 30, shade = .1, border = NA)
#par(new=TRUE)
persp(bivn.kde2, phi = 45, theta = 30, shade = .1, border = NA)
```

<br><br>

#### <font color = "#880000">Comment</font>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
Remmeber, that once we know the joint distribution function all marginals are easily obtainable. On the other hand, however, once we know all marginals it is still not enough to reconstruct the whole joint (multivariate) distribution. Using this idea we can do some the next principle to generate some multivariate normal distribution BUT <font color ="red">we can not control for the overal dependence structure - the covariance matrix among the marginals.</font>

```{r, message=F, warning=FALSE}
library(ggplot2)
```

```{r}
rm(list = ls())
set.seed(1234)

n <- 1000
sample1 <- rnorm(n, mean=2, sd = 2)
sample2 <- 0.2 * sample1 + rnorm(n)
data <- data.frame(x = sample1, y = sample2, group="A")

sample3 <- rnorm(n, mean=0, sd = 1)
sample4 <- 3*sample3 - 1 + rnorm(n)
data <- rbind(data, data.frame(x = sample3, y = sample4, group="B"))
```

```{r, fig.width =10}
ggplot(data, aes(x=x, y=y,colour=group)) + geom_point(size=1.5, alpha=.6) 
```


<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: third lab session / 21.03.2016)

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 contain a scatterplot 
of the two covariates (see the previous item) including the contours which you calculate from the sample mean vector and the sample variance-covariance matrix.</li>

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

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