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

## <b>Conditional and marginal distributions of random vectors</b>
### <font color = "#880000"><b>Summer Term 2021/2022 |  22/02/22 </b></font>
###### <a href = "Lab2_2022.Rmd">sourse file (UTF8 encoding)</a>

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

Outline of the second NMST539 lab session:

<ul>
<li>matrix eigen decomposition, eigen values, and eigen vectors;</li>
<li>random vectors: conditional and marginal distributions;</li>
<li>multivariate normal distribution with R packages;</li>
</ul>

The first task of the lab session today is to recall some important matrix properties (e.g., eigen value decomposition, eigen values properties, eigen vectors) crutial 
for most of the multivariate statistical techniques discussed later (during the term). Some teoretical/practical illustrations will be used and tge multivariate normal distribution
and its specific implementation in R will be discussed as well. 

The R-software is available for download (GNU licence, free of charge) from the website: <a href = "https://www.r-project.org">https://www.r-project.org</a>

A user-friendly interface (one of many) however, not nesessary for a proper functioning of the R software: <a href = "https://www.rstudio.com">RStudio</a>.<br><br>

<b>Some (rather brief) manuals and introduction tutorials for <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>

<b>Theoretical bibliography:</b><br>
<ul>
<li>Hardle, W. and Simar.L.: Applied Multivariate Statistical Analysis. Springer, 2015</li>
<li>Mardia, K., Kent, J., and Bibby, J.:Multivariate Analysis, Academic Press, 1979.</li>
</ul><br>

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



### <font color = "#880000"><b>1. Eigen values and eigen vectors</b></font>
Multivariate statistical methods are mostly based on some proper matrix decomposition. Usually we are either dealing with a decomposition of a rectangular data matrix -- a data frame 
of a type $n \times p$ or, instead, a square type matrix ($p \times p$) such as, for instnace, a theoretical variance-covariance matrix, its 
sample couterpart respectively, or the precision/distance matrix.

The most important decomposition is the <b>eigen value matrix decomposition</b> (i.e., eigendecomposition or a spectral decomposition) of a square type matrix into a set of so called 
eigen vectors and the correponding eigen values. 

Nonzero vector $\boldsymbol{v} \in \mathbb{R}^p$ is said to be an eigen vector of the matrix $\boldsymbol{A} \in \mathbb{R}^{p \times p}$ if 
$$
\boldsymbol{A}\boldsymbol{v} = \lambda \boldsymbol{v},
$$
for   $\lambda \in \mathbb{R}$ being some fixed number -- the<b> eigen value of the matrix $\boldsymbol{A}$</b> corresponding to the <b>eigen vector $\boldsymbol{v} \in \mathbb{R}$</b>. 

<ul>
<li>Interesting illustration of the role of the eigen vectors and eigen values can be fond here:<br>
<a href = "http://setosa.io/ev/eigenvectors-and-eigenvalues/">http://setosa.io/ev/eigenvectors-and-eigenvalues/</a></li>

<li>Some rather brief and basic information can be also found, for instance, here:<br> <a href = "https://personal.utdallas.edu/~Herve/Abdi-EVD2007-pretty.pdf">https://personal.utdallas.edu/~Herve/Abdi-EVD2007-pretty.pdf</a></li>
</ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />



#### <font color = "#880000">Individual task</font>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<ul>
<li>Recall some important matrix properties and the corresponding facts about the eigen values and eigen vectors. </li>

<li>Try to make your own R code which will be able to find the corresponding eigen values of a simple 
$2 \times 2$ or $3 \times 3$ square (symmetric) matrix. </li>

<li>The spectral decomposition of a square matrix (of a type $p \times p$) in R (eigendecompostion - EIV) is implemented within a function `eigen()`.
Use the help and learn how the function is supposed to be applied to find EIV decomposition of some real (square) matrix. 
</li>

<li>The EIV decomposition of a square matrix  $\boldsymbol{A} \in \mathbb{R}^{p \times p}$ is generally defined as 
$$
 \boldsymbol{A} = \boldsymbol{Q}\boldsymbol{\Lambda}\boldsymbol{Q}^{-1}
$$
where $\boldsymbol{\Lambda} \in \mathbb{R}^{p \times p}$ is a diagonal matrix with the corresponding eigen values on its diagonal. 
Matrix $\boldsymbol{Q} \in \mathbb{R}^{p \times p}$ contains the eigen vectors (as columns). </li>
</ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />

<br><br>



### <font color = "#880000"><b>2. Random vectors (join, marginal, conditional distributions)</b></font>

Multivariate statistical methods are based on the stochastic theory of random vectors (random matrices respecticely). 
In this part we use a few illustrative (theoretical and practical) examples to highlight some crucial properties. The examples will 
be also used to illustrate some important approaches used for simulating multivariate random samples (a general algorithm applicable not only for the R programming). 


<b>Theoretical and practical examples</b>
<ol>
<li>Consider a uniform (two dimensional) distribution of a random vector $(X,Y)^\top$ over the set $M = \{(x, y) \in \mathbb{R}^2;~0 < x < y < 1\}$. </li>
Find the explicit expression for the joint distribution -- the density function $f_{(X,Y)}(x,y)$ of the random vector $(X, Y)^\top$. Derive also the corresponding 
marginal densities $f_X(x)$ and $f_Y(y)$. Can you say something relevant about the dependence/independence of the random variables  $X$
and $Y$ (elements of the random vector $(X, Y)^\top$)? </li><br>

<li>Use the expressions above to simulate a two dimensional random sample $(X_1, Y_1)^\top, \dots, (X_n, Y_n)^\top$ 
from the joint distribution given by the two-dimensional density function $f_{(X,Y)}(x,y)$ (from the previous example). 
Use some proper graphical tools to verify that the simulated random sample indeed follows the given joint distribution. 
Calculate important sample characteristics (sample means  and the sample variance-covariance matrix). <br><br>

<b>Practical hint</b>
<ul>
<li>For an arbitrary random variable $X$ with the distribution function $F_{X}(x)$ we have that 
$$
P[X \leq x] = F_{X}(x),
$$
for all $x \in \mathbb{R}$. If, in addtion, the distribution function is absolutely continuous, then for the quantile function $F_X^{-1}(x)$ it holds that it is 
obtained as an inverse function of the distribution function  $F_{X}(x)$. Thus, for the random variable  $Z = F_X(X)$ it holds  (following the definition of the distribution function and quantile function), that
$$
P[Z \leq z] = P[F_X(X) \leq z] = P[X \leq F_X^{-1}(z)] = F_{X}\Big( F_X^{-1}(z)\Big) = z
$$
for $z$ from the unit interval $[0,1]$. This means, that the random variable $Z$ follows a uniform distribution on the interval $[0,1]$. <br>

This should be enough to generate a random sample basically from an arbitrary distribution (if the distribution function is known). <br>An illustrative example 
is given below. 
</li><br>

<li>Consider a random variable $X$ with the distribution given by the uniform density function $f(x) = 3 x^2$ for $x \in (0,1)$ and $f(x) = 0$ otherwise.
We want to use the R program to simulate a random sample from such distribution -- however, without recalling a fact that actually the random variable  $X$ follows the Beta distribution
with parameters $\alpha = 3$ and $\beta = 1$. 

The distribution function of the random variable $X$ is $F(x) = x^3$
for $x \in [0,1]$ (while the distribution function is zero for  $x < 0$ and $F(x) = 1$ for $x > 1$). As far as the distribution function is absolutely continuous
we can set the corresponding quantile function  $F^{-1}(q)$, for $q \in [0,1]$ as the inverse function of $F(x)$, i.e., 
$$
F^{-1}(q) = \sqrt[3]{q}.
$$
for $q \in [0,1]$. Of course, it holds that $F^{-1}(0) = 0$ and  $F^{-1}(1) = 1$. Thus, the random variable $Z = F(X)$, for $F(x) = x^3$ being the distribution function of $X$,
follows a uniform distribution on the interval $[0,1]$. Therefore, we can use a simple (pseudo)-random generator in R (function `runif()`) to generate a random sample $Z_1, \dots, Z_n$ from the uniform distribution on the interval $[0,1]$.

```{r, fig.width = 10}
n <- 1000
set.seed(1234)

Z <- runif(n)
par(mfrow = c(1,2))
plot(ecdf(Z), main = "Empirical d.f. of Z")
hist(Z, col = "gray", breaks = ceiling(sqrt(n)))
```

Next, we will use the back-transformation $X_i = F^{-1}(Z_i)$ in order to obtain the random variables $X_i$ which will already follow the distribution given by the distribution function $F_X(x)$, or the density function $f(x)$ respectively. 

```{r, fig.width = 10}
X <- Z^(1/3)
xSeq <- seq(0,1, length = 1000)
par(mfrow = c(1,2))
plot(ecdf(X), main = "Empirical d.f. of X")
lines(xSeq^3 ~ xSeq, col = "red", lty = 2, lwd = 2)
legend("topleft", legend = c("Empirical d.f.", "Theoretical d.f."), col = c("black", "red"), lty = c(1,2), lwd = c(2,2))
hist(X, freq = F, col = "gray", breaks = ceiling(sqrt(n)))
lines(3 * xSeq^2 ~ xSeq, col = "red", lty = 2, lwd = 2)
legend("topleft", legend = c("Empirical density", "Theoretical density"), col = c("black", "red"), lty = c(1,2), lwd = c(2,2))
```
</li>

<li>Use the previous example for a larger sample size $n \in \mathbb{N}$ and use some empirical tools to verify that the given approach 
indeed provides a sample form the right distribution. The number of bins in the histogram figure can be set by an additional parameter  `breaks = ...`. </li>

<li>In order to use the approach described above for generating a two dimensional random sample $(X_i,Y_i)^{\top}$ for $i = 1, \dots, n$ (from some given joint distribution)
it is enough to apply analogous steps separatelly  for the random element  $X$ and also for $Y$ (if these two are stochastically independent) or, instead, it is 
needed to use the analogous steps for the first element only and, later, to apply analogous steps considering the conditional distribution of the second element given the first (if the variables are stochastically dependent).</li>
</ul></li><br>

<li>
Use the two-dimensional Pareto distribution given by the joint density function  
$$
f(x, y) = c(x + y - 1)^{-p - 2}, 
$$
for $x,y > 1$ and $p > 0$. Find the value of the normalizing constant  $c > 0$ and the explicit formulas for the marginal distributions of both random elements. 
Derive the conditional distribution of 
$X$ given $Y = y$ and the conditional distribution of $Y$ given $X = x$. 
</li>
</ol>

<br><br>

#### <font color = "#880000">Individual task</font>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<ul>
<li>How would you use the approach described above to visually verify independence/dependence of some given (two) random variables from the finite 
 sample $(X_i, Y_i)^\top$, for $i = 1, \dots, n$.</li>

<li>A standard  $xy$ scatterplot is usually the most common choise, howver, very often it is  not an appropriate tool to judge about dependence/independence of two random variables.  See, for instance, 
the figure below: 

```{r, fig.width=10}
n <- 100

set.seed(1234)
x <- rexp(n, rate = 1)
y <- rexp(n, rate = 10)
z <- x + rnorm(n, -1, 2)

par(mfrow = c(1,2))
plot(y ~ x, pch = 21, bg = "lightblue")
plot(z ~ x, pch = 21, bg = "lightblue")
```

The first figure shows two independent random variables $X$ and $Y$. Despite the fact that there are some clear clusters of points near the zero point $[0,0]$, the random variables $X$ and $Y$
are indeed stochastically independent. 

On the other hand, the second figure shows the scatterplot for two dependent random variables    $X$ and $Z$ -- however, using just the scatterplot itself, maybe we would have a tendency two consider $X$ and $Z$ to be rather independent. 

The whole picture, however, changes when we consider different scatterplots -- instead of plotting the original data $X_i$, $Y_i$, and $Z_i$ (for $i = 1, \dots, n$) 
we can plot analogous scatterplots for $\widehat{F}_{X}(X_i)$, $\widehat{F}_{Y}(Y_i)$, and 
$\widehat{F}_{Z}(Z_i)$, where $F_{X}$, $F_{Y}$, and $F_{Z}$ are the corresponding distribution functions of $X$, $Y$, and $Z$.


```{r, fig.width=10}
qx <- NULL 
qy <- NULL 
qz <- NULL 

for (i in 1:n){
   qx <- c(qx, which(sort(x) == x[i]))
   qy <- c(qy, which(sort(y) == y[i]))
   qz <- c(qz, which(sort(z) == z[i]))
} 

qx <- qx / n
qy <- qy / n
qz <- qz / n

par(mfrow = c(1,2))
plot(qy ~ qx, pch = 21, bg = "lightblue", xlab = "", ylab = "")
lines(lowess(qy ~ qx, f = 1), col = "red", lwd = 2)
plot(qz ~ qx, pch = 21, bg = "lightblue", xlab = "", ylab = "")
lines(lowess(qz ~ qx, f = 1), col = "red", lwd = 2)
```

</li>

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

<br><br>





### <font color = "#880000"><b>3. Multivariate normal distribution</b></font>

<font color = "#880000"><b>Library `mvtnorm`</b></font>

Firstly, we need to make sure that the library for the multivariate normal distribution (R package `mvtnorm`) is installed and properly 
loaded in the R working environment (use commad `library("mvtnorm")`). If the library can not be loaded you need to install it first. The installation is 
performed fully automatically by typing the command `install.packages("mvtnorm")` in the R command line. 

```{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$ (a symmetric and positive definite matrix). 

For more details about the library (including a PDF 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 some specific purposes (master thesis, scientific paper, simulations, etc.) this is even a strictly required step. 

For instance try (among each other) the command:

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

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

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

Is it clear, what is the difference between these two outputs? <br>
Producing results which are easily reconstructable and can be verified is the key issue here. 

<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">Example</font>

Let us start with a simple (two-dimensional) example: we are going to generate a random sample of size $n \in \mathbb{N}$ from the two-dimensional normal distribution

<center>
$\boldsymbol{X}_{1}, \dots, \boldsymbol{X}_n \sim 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)$.
</center>
for $\boldsymbol{X}_{i} = (X_{i 1}, X_{i 2})^\top$.

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

<ul>
<li>What is the marginal distribution of the elements -- random variables $X_{i 1}$ and $X_{i 2}$?</li>
<li>Can we use the standard -- one dimensional random generator `rnorm()` to simulate both variables individually and, later, to combine them into the vector 
to form the multivariate -- two dimensional normal distribution?
<li>Verify by yourself that the variance-covariance matrix $\Sigma$ is a symmetric and positive definite matrix.</li>
</ul>

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 standard R command `plot()`:

```{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", cex = 0.8)
```

Using now the 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 expression for 
the mutivariate normal density, 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 generated 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 defined as

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

where $\boldsymbol{x} = (x_{1}, x_{2})^\top \in \mathbb{R}^2$, $\boldsymbol{\mu} = (\mu_{1}, \mu_{2})^\top$ is the mean vector, and $\Sigma$ is some given variance-covariance matrix (which is symmetric and positive definite).

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 now applied in the following way:

```{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", cex = 0.8)
contour(x,y,matrix(z,n,n), levels=(0:15), col = terrain.colors(16), add=TRUE, lwd = 1)
```


<br><br>

#### <font color = "#880000">Individual task</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 values of the given 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><br>

<li>Can you think of some reasonable and meaningful interpretation of the plotted contours in the figure above? Why are they called <i>depth contours</i>? Is there some straightforward 
relationship with quantiles?</li><br>

<li>Think of some linear combination of the random elements $X_{i1}$ and $X_{i2}$ in the sample you generated from the two dimensional normal distribution. What is the distribution of this linear combination? Use some graphical tools to verify your conclussion. </li><br>
</ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />

<br><br>


<font color = "#880000"><b>Library `mixtools`</b></font>

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">Individual task</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>Library `MASS`</b></font>

One more 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 find out more details about this command). 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, for some visualisation purposes, we calculate the corresponding 
multivariate kernel density estimates (theoretical details are omitted and they are not important). 

```{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 -- one just need to calculate some multivariate integrals which is usually easily doable. On the other hand, however, if we know all marginals it is still not enough to reconstruct the whole joint (multivariate) distribution. Using this idea we can use 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) 
```


Checking independence by using the distribution functions: 

```{r, fig.width =10, fig.height=10}
xPoints <- pnorm(sample1, mean(sample1), sd(sample1))
yPoints <- pnorm(sample2, mean(sample2), sd(sample2))

plot(xPoints, yPoints, pch = 21, bg = "lightblue")
lines(lowess(yPoints ~ xPoints), col = "red", lwd = 2)
```



Which can be compared with the figure for plotting sample no.1 against the sample no.3 (two independent samples). 

```{r, fig.width =10, fig.height = 10}
xPoints <- pnorm(sample1, mean(sample1), sd(sample1))
yPoints <- pnorm(sample3, mean(sample3), sd(sample3))

plot(xPoints, yPoints, pch = 21, bg = "lightblue")
lines(lowess(yPoints ~ xPoints), col = "red", lwd = 2)
```


#### <font color = "#880000">Comments & Questions</font>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<ul>
<li>What do you expect to see in the two figures displayed above?</li>

<li>Do you understand the idea of ploting the values of the distribution functions evaluated at the random samples against each other?</li>

<li>Which distributions can be used to this approach? </li>

<li>Why is it not a good idea to plot directly the values of the sample against each other?</li>

</ul>
 







<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>Individual task assignment</b></font>
### (Deadline: Third lab session / 01.03.2022)

<ul>
<li> In the second section ''Random vectors (join, marginal, conditional distributions)''  solve Example 1 (theoretical) 
and use the R software to simulate the random sample in Example 2 (practical).</li><br>

<li> Place the solution of both examples on you website, on<font color = "red">Tuesday, 01.03.2022, at 12:00 at lastest</font>.</li>
</ul>

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