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

## <b>Wishart and Hotelling's Distribution</b>
###  (application for confidence bands and statistical tests)
### <b><font color = "#880000">Summer Term 2018/2019 | Monday 11/03/19</font></b>
###### <a href = "cvicenie2019_4.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. Wishart's distribution</b></font>

The Wishart distribution is a multivariate generalization of the univariate $\chi^2$ distribution used for the variance inference in the univariate case. It is named in honor of John Wishart, who  formulated this distribution in 1928. The Wishart distribution is used as an eficient tool for the analysis of the variance-covariance matrix of same random sample $\boldsymbol{X}_{1}, \dots, \boldsymbol{X}_{n}$, for some $n \in \mathbb{N}$, where  $\boldsymbol{X}_{i} = (X_{i 1}, \dots, X_{i p})^\top$ is a $p$-dimensional random vector (for instance, $p$ different covariates/characteristics recorded on each subject). 



The Wishart distribution is a multi-variate generalization of the $\chi^2$ distribution in the following sence: for instance, for a normally distributed random sample $\boldsymbol{X}_{1}, \dots, \boldsymbol{X}_{n}$ drawn from some multivariate distribution $N_{p}(\boldsymbol{0}, \Sigma)$, with the zero mean vector and $\Sigma$ being some variance-covariance matrix (symmetric and positive definite), we have that 
$$
\mathbb{X}^{\top}\mathbb{X} \sim W_{p}(\Sigma, n), 
$$
for $\mathbb{X} = (\boldsymbol{X}_{1}^\top, \dots, \boldsymbol{X}_{n}^\top)^{\top}$. Similarly, for a univariate random sample $X_{1}, \dots, X_{n}$ drawn, for instance, from $N(0,1)$, we have that 
$$
\boldsymbol{X}^{\top}\boldsymbol{X} \sim \chi_{n}^2, 
$$
where $\boldsymbol{X} = (X_{1}, \dots, X_{n})^\top$. Thus, for a sample of size $n \in \mathbb{N}$ drawn from $N(0,1)$ the corresponding distribution of $\mathbb{X}^\top\mathbb{X}$  is $W_{1}(1, n) \equiv \chi_{n}^2$.  Thus, the Wishart distribution is a family of probability distributions defined over symmetric, nonnegative-definite matrix-valued random variables (so called <i>random matrices</i>). 

The corresponding density function of the Wishart distribution takes the form
$$
f(\mathcal{X}) = \frac{1}{2^{np/2} |\Sigma|^{n/2} \Gamma_{p}(\frac{n}{2})} \cdot |\mathcal{X}|^{\frac{n - p - 1}{2}} e^{-(1/2) tr(\Sigma^{-1}\mathcal{X})}, 
$$
for $\mathcal{X}$ being a $p\times p$ random matrix and $\Gamma_{p}(\cdot)$ is a multivariate generalization of the Gamma function $\Gamma(\cdot)$
In the R software there are various options (packages and commands) on how to use and apply the Wishart distribution. 

<br><br>

#### <font color = "#880000"><b>To Do by Yourself (Theoretical and Practical)</b></font>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<ul>
<li>Consider a  multivariate random sample $\boldsymbol{X}_{1}, \dots, \boldsymbol{X}_{n}$ which comes from some general multivariate distribution $N_{p}(\boldsymbol{\mu}, \Sigma)$, with some mean vector $\boldsymbol{\mu} \in \mathbb{R}^p$ and the variance-covariance positive definite matrix $\Sigma$. Show, that the random matrix defined as $n\mathbb{X}^\top \mathcal{H}\mathbb{X}$ follows the Wishart distribution $W_{p}(\Sigma, n - 1)$, for $\mathcal{H}$ being the idenpotent centering matrix $\mathcal{H} = \mathbb{I}_{n} - \frac{1}{n}\boldsymbol{1}_{n}\boldsymbol{1}_{n}^\top$.</li><br>

<li>use the help session in R to find out more about some standard commands used in R for the Wishart distribution (in particular, check the help session for the command `rWishart()`);</li><br>

<li>some other commands and additional extensions are available  in some packages which need to be downloaded and installed separately, e.g. function `Wishart()` in the library 'MCMCpack', or `Wishart()` in the library 'mixAK', or `rwishart()` in the library 'dlm', `dist.Wishart()` in the library 'LaplacesDemon', and others;</li>
</ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<br>

A simple random sample from the univariate Wishart distribution $W_{1}(\Sigma = 1, n = 10)$ (equivalently a $\chi^2$ distribution with $n = 10$ degrees of freedom) can be obtained, for instance, as 

```{r}
S <- as.matrix(1)
sample <- rWishart(10, df = 10, S)
```

For a general $p \in \mathbb{N}$ the same command can be used, however, with the corresponding variance-covariance matrix $S$ to be properly specified 
as a symetric, positive nefinite matrix $\Sigma$, for instance:

```{r}
Sigma <- cbind(c(1,0,0,0), c(0,1,0,0), c(0,0,1,0), c(0,0,0,1))
rWishart(2, df = 100, Sigma) ### sample of size two
```


We can also use the standard approach for generating a sample from the $\chi^2$ distribution -- the R function `rchisq()` --  and, using some histograms, we can empirically  compare 
the univariate Wishart distribution and the $\chi^2$ distribution. For large sample size $n \in \mathbb{N}$ we they should be quite similar/identical. 

```{r, fig.width = 10}
set.seed(1234)
sampleWishart <- rWishart(5000, df = 10, S)
sampleChiSq <- rchisq(5000, df = 10)

par(mfrow = c(1,2))
hist(sampleWishart, col = "lightblue", main = expression(paste("Wishart Distribution", sep ="")), xlim = c(0, 40), breaks = 30, freq = F)
lines(density(sampleWishart), col = "red", lwd = 2)
hist(sampleChiSq, col = "lightblue", main = expression(paste(chi^2, "Distribution", sep = "")), xlim = c(0,40), breaks = 30, freq = F)
lines(density(sampleChiSq), col = "red", lwd = 2)
```


<br><br>

#### <font color = "#880000"><b>To Do by Yourself</b></font>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<ul>
<li>Consider some Wishart distribution for some $p \in \mathbb{N}$ such that $p > 1$ and generate a random sample from the Wishart distribution $W_{p}$ 
with some reasonable parameters $\Sigma$ and $df$.</li><br>

<li>Consider some general matrix $A \in \mathbb{R}^{p \times q}$ for some $q \in \mathbb{N}$ such that $p \neq q$ and use some visualization tools to empiricaly verify (for instance, by some histograms, or KDE) that $A^\top \mathbb{M} A \sim W_{q}(A^\top\Sigma A, n)$, where $\mathbb{M} \sim W_{p}(\Sigma, n)$; </li>
</ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<br><br>

#### <font color = "#880000"><b>2. Hotelling's $\boldsymbol{T^2}$ Distribution</b></font>

In a similar manner as we define a classical $t$-distribution in a univariate case (i.e. standard normal $N(0,1)$ variable devided by a square root of a $\chi^2$ variable normalized by its degrees of freedom) we define a multivariate generalization (the Hotelling's $T^{2}$ distribution) as
$$
n \boldsymbol{Y}^{\top} \mathbb{M}^{-1} \boldsymbol{Y} \sim T^{2}(n, p),
$$
where $p \in \mathbb{N}$ is the dimension of the random vector $Y \sim N_{p}(0, \mathbb{I})$ and $n \in \mathbb{N}$ being the parameter of the Wishart distribution for the random matrix $\mathbb{M} \sim W_{p}(\mathbb{I}, n)$. Moreover, the random vector $\boldsymbol{Y}$ is assumed to be independent of the random matrix $\mathbb{M}$.



A special case for $p = 1$ gives the standard Fisher distribution with one and $n$ degrees of freedom (equivalently a square of the $t$-distribution with $n$ degrees of freedom). However, the Hotelling's $T^2$ distribution with parameters $p, n \in \mathbb{N}$ is, in general, closely related with  the Fisher's F distribution. It holds, that 

$$
T^{2}(p, n) \equiv \frac{n p}{n - p + 1}F_{p, n - p + 1}.
$$

Therefore, the standard univariate Fisher's F distribution can be effectively used to draw critical values also from the Hotelling's $T^2$ distribution. 
The corresponding transformation between the Hotelling's $T^2$ distribution and the Fisher's F distribution only depends on parameters $n, p \in \mathbb{N}$. <br><br>


The role of the given parameters and different parameter settings in the Hotelling's $T^2$ distribution can be (a little) visualized, for instance, by utilizing the Fisher's F distribution
and the knowledge about its mean and variance: 

```{r, fig.width = 10, fig.height=10}
samples <- NULL
samples <- cbind(samples, rf(n = 1000, df1 = 1, df2 = 1))
samples <- cbind(samples, rf(n = 1000, df1 = 1, df2 = 10))

samples <- cbind(samples, rf(n = 1000, df1 = 10, df2 = 10))
samples <- cbind(samples, rf(n = 1000, df1 = 10, df2 = 100))

samples <- cbind(samples, rf(n = 1000, df1 = 100, df2 = 10))
samples <- cbind(samples, rf(n = 1000, df1 = 100, df2 = 100))

samples <- cbind(samples, rf(n = 1000, df1 = 1000, df2 = 100))
samples <- cbind(samples, rf(n = 1000, df1 = 1000, df2 = 1000))

par(mfrow = c(4,2))
for (i in 1:dim(samples)[2]){
  hist(samples[,i], xlab=paste("Sample no.", i, sep = ""), col = "yellow", main = "", freq = F, breaks = 20)
  lines(density(samples[,i]), col = "red")
}
```

Let us recall, that for some random variable $X$ with the standard Fisher's F distribution $F_{df_1, df_2}$ we have the following: 

<ul>
<li> expected value $E X = \frac{df_{2}}{df_{2} - 2}$, for $df_2 > 2$;</li>
<li> variance $Var(X) = \frac{2 df_{2}^2 (df_1 + df_2 - 2)}{df_1 (df_{2} - 2)^2 (df_2 - 4)}$ for $df_2> 4$;</li>
</ul>

<br><br>
Moreover, using the relationship between the Hotelling's $T^2$ distribution and the Fisher's F distriubtion we can effectively use the $F$-distribution to draw the corresponding critical values when needed for some statistical tests. 

For some multivariate random sample $\boldsymbol{X}_{1}, \dots, \boldsymbol{X}_{n} \sim N_{p}(\boldsymbol{\mu}, \Sigma)$ for some unknown mean vector $\boldsymbol{\mu} \in \mathbb{R}^{p}$ and some variance-covariance matrix $\Sigma$, we have that 
$$
(n - 1)\Big(\overline{\boldsymbol{X}} - \boldsymbol{\mu}\Big)^\top \mathcal{S}^{-1} \Big(\overline{\boldsymbol{X}} - \boldsymbol{\mu}\Big) \sim T^2(p, n - 1), 
$$
which can be also equivalently expressed as 
$$
\frac{n - p}{p} \Big(\overline{\boldsymbol{X}} - \boldsymbol{\mu}\Big)^\top \mathcal{S}^{-1} \Big(\overline{\boldsymbol{X}} - \boldsymbol{\mu}\Big) \sim F_{p, n - p}. 
$$
This can be now used to construct confidence regions for the unknown mean vector $\boldsymbol{\mu}$ for testing hypothesis about the true value of the vector of parameters $\boldsymbol{\mu} = (\mu_{1}, \dots, \mu_{p})^{\top}$. However, rather than constructing a confidence region for $\boldsymbol{\mu}$ (which can be impractical in higher dimensions for even slightly larger values of $p$) one focusses on construction confidence intervals for the elements of $\boldsymbol{\mu}$ such that the mutual coverage is under control (usually we require a simultaneous coverage of $(1 - \alpha)\times 100~\%$ for some small $\alpha \in (0,1)$).

For a hypothesis test 

$$
H_{0}: \boldsymbol{\mu} = \boldsymbol{\mu}_{0} \in \mathbb{R}^{p}
$$

$$
H_{1}: \boldsymbol{\mu} \neq \boldsymbol{\mu}_{0} \in \mathbb{R}^{p}
$$

we can use the following test statistic:

$$
(n - 1)\Big(\overline{\boldsymbol{X}} - \boldsymbol{\mu}_{0}\Big)^\top \mathcal{S}^{-1} \Big(\overline{\boldsymbol{X}} - \boldsymbol{\mu}_{0}\Big), 
$$

which, under the null hypothesis, follows the $T^2(p, n - 1)$ distribution. Equivalently we also have that 

$$
\frac{n - p}{p} \Big(\overline{\boldsymbol{X}} - \boldsymbol{\mu}_{0}\Big)^\top \mathcal{S}^{-1} \Big(\overline{\boldsymbol{X}} - \boldsymbol{\mu}_{0}\Big) 
$$

follows the Fisher $F$-distribution with $p$ and $n - p$ degrees of freedom.


In the R software we can use the library `DescTools` (the library needs to be installed by running the command `install.packages("DescTools")`) which is loaded into  the R environment by the command 

```{r}
library(DescTools)
```

Now we can use the function `HotellingsT2Test()` to perform the test on multivariate mean vector (use the help session to find out how the function works). 


Using the same approach we can also construct the confidence elipsoid for $\boldsymbol{\mu} \in \mathbb{R}^p$ - it holds that 
$$
\frac{n - p}{p} \Big(\overline{\boldsymbol{X}} - \boldsymbol{\mu}_{0}\Big)^\top \mathcal{S}^{-1} \Big(\overline{\boldsymbol{X}} - \boldsymbol{\mu}_{0}\Big)  \sim F_{p, n - p},
$$
and therefore, the following set 
$$
\left\{\boldsymbol{\mu} \in \mathbb{R}^p;~  \Big( \overline{\boldsymbol{X}} - \boldsymbol{\mu}\Big)^\top \mathcal{S}^{-1} \Big( \overline{\boldsymbol{X}} - \boldsymbol{\mu}\Big) \leq \frac{p}{n - p} F_{p, n- p}(1 - \alpha) \right\}
$$
is a confidence region at the confidence level of $\alpha = 1 - \alpha$ for the vector of parameters $\boldsymbol{\mu} \in \mathbb{R}^p$ - an interior of iso-distance ellipsoid in $\mathbb{R}^p$.

A brief example on how it works (example from the lecture notes): 
```{r, fig.width=10}
library(mvtnorm)
s=matrix(c(1,-0.5,-0.5,1),2);x=seq(-3,3,by=0.015) ### variance-covariance matrix

contour(x,x,outer(x,x, function(x,y){dmvnorm(cbind(x,y),sigma=s)}))
n <- 20
X <- rmvnorm(n,sigma=s)
m <- apply(X,2,mean)
S <- cov(X)
points(m[1],m[2],pch=8,col="red",cex=2)

S1=solve(S)
contour(x,x,outer(x,x,function(x,y){(n-2)*
                                      apply(t(t(cbind(x,y))-m),1,function(x){t(x)%*%S1%*%x})<
                                      2*qf(0.95,2,n-2)}),col="red",add=TRUE)
bodyx=m[1]+c(-1,1)*sqrt(S[1,1]*2*qf (0.95,2,n-2) /(n-2))
bodyy=m[2]+c(-1,1)*sqrt(S[2,2]*2*qf (0.95,2,n-2) /(n-2))
polygon(x=bodyx[c(1,1,2,2,1)],y=bodyy[c(1,2,2,1,1)],border="blue")
```

(which is the 95% confidence region the true mean vector)<br><br>


#### <font color = "#880000"><b>To Do by Yourself</b></font>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<ul>
<li>Use the approach described above and construct a confidence region for the true mean vector for some data set of your choice;</li>
<li>Explain and interpret the obtained results. What is the conclusion?</li>
</ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<br><br>



#### <font color = "#880000"><b>Difference of two multivariate means with the same variance-covariance matrix</b></font>
In an analogous way one can also construct a two sample Hotelling test to compare two population means. 

We assume a multivariate sample $\boldsymbol{X}_{1}, \dots, \boldsymbol{X}_{n} \sim N_{p}(\boldsymbol{\mu}_{1}, \Sigma)$
and some another sample $\boldsymbol{Y}_{1}, \dots, \boldsymbol{Y}_{M} \sim N_{p}(\boldsymbol{\mu}_{2}, \Sigma)$, with some generally different mean parameter $\boldsymbol{\mu}_{1} \neq \boldsymbol{\mu}_{2}$. We are interested 
in testing the null hypothesis 
$$
H_{0}: \boldsymbol{\mu} = \boldsymbol{\mu}_{0}
$$
against the alternative hypothesis, that the null does not hold. We have that the following holds: 
$$
(\overline{X}_{1} - \overline{X}_{2}) \sim N_{p}\Bigg(\boldsymbol{\Delta}, \frac{n + m}{n m} \Sigma\Bigg),
 $$
and also 
$$
n\mathcal{S}_{1} + m\mathcal{S}_{2} \sim W_{p}(\Sigma, n + m - 2), 
$$
where $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$ respectively are the empirical estimates of the variance-covariance matrix $\Sigma$, given the first sample and the second sample respectively. Then the rejection region is defined as 
$$
\frac{nm(n + m - p - 1)}{p(n + m)^2}(\overline{\boldsymbol{X}}_{1} - \overline{\boldsymbol{X}}_{2})^{\top} \mathcal{S}^{-1} (\overline{\boldsymbol{X}}_{1} - \overline{\boldsymbol{X}}_{2}) \geq F_{p, n + m - p - 1}(1 - \alpha).
$$
<br>
<hr size="3" width="100%" noshade style="color:#440000" align="left" /><br><br>

#### <font color = "#880000"><b>Difference of two multivariate means with unequal variance-covariance matrix</b></font>
Now, we assume a multivariate sample $\boldsymbol{X}_{1}, \dots, \boldsymbol{X}_{n} \sim N_{p}(\boldsymbol{\mu}_{1}, \Sigma_{1})$
and some another sample $\boldsymbol{Y}_{1}, \dots, \boldsymbol{Y}_{M} \sim N_{p}(\boldsymbol{\mu}_{2}, \Sigma_{2})$, with some generally different mean parameter $\boldsymbol{\mu}_{1} \neq \boldsymbol{\mu}_{2}$. We are interested 
in testing the null hypothesis 
$$
H_{0}: \boldsymbol{\mu} = \boldsymbol{\mu}_{0}
$$
against the alternative hypothesis, that the null does not hold. Again, we have that the following holds: 
$$
(\overline{X}_{1} - \overline{X}_{2}) \sim N_{p}\Bigg(\boldsymbol{\Delta}, \frac{\Sigma_{1}}{n} + \frac{\Sigma_{2}}{m}\Bigg),
 $$
and therefore, it also holds that 
$$
(\overline{X}_{1} - \overline{X}_{2})^\top \Big(\frac{\Sigma_{1}}{n} + \frac{\Sigma_{2}}{m}\Big)^{-1} (\overline{X}_{1} - \overline{X}_{2}) \sim \chi_{p}^{2}.
$$


<br><br>

#### <font color = "#880000"><b>To Do by Yourself</b></font>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<ul>
<li>download and install the R library 'Hotelling' and 'DescTools';</li>
<li>for one sample and two sample Hotelling tests use the command `HotellingsT2Test()` from 'DescTools' package;</li>
<li>for two sample Hotelling tests you can also use the command `hotelling.test()` from 'Hotelling' package;</li>
<li>use the help session to find out more about commands `hotelling.stat()` and `hotelling.test()`
which can be used to test the null hypothesis about the true value of the mean vector;
</ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<br><br>

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

There is a dataset called `container.df` available in the package 'Hotelling'. The data represent some concentration measurements on nine different elements (see the help session `?container.df` for further details). 

```{r}
data(container.df)
summary(container.df)
```

The estimates for the mean vector and variance-covariance matrix can be obtained as 
```{r}
apply(container.df[,2:10], 2, mean) ### mean vector estimate
cov(container.df[,2:10])            ### variance-covariance estimate
```
and we can test whether the mean vector of concentrations in the first container (`gp == 1`) and 
the second container (`gp == 2`) is the same or not. Try to perform a one sample Hotelling test to 
test whether the unknown vector of mean parameters equals some prespecified vector $\boldsymbol{\mu}_{0} \in \mathbb{R}^{p}$.<br><br>



For a two sample problem we can use the first covariate in the data to define two populations and we provide a test whether the mean vectors are equal or not. The corresponding Hotelling test statistic equals

```{r}
split.data = split(container.df[,-1],container.df$gp)
x <- split.data[[1]]
y <- split.data[[2]]
statistic <- hotelling.stat(x, y)
```

and the corresponding test is performed by

```{r}
hotellingTest <- hotelling.test(.~gp, data = container.df)
hotellingTest
```

<br><br>

#### <font color = "#880000"><b>Questions</b></font>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<ul>
<li>What is the interpretation of the test results?</li>
<li>Try to perform a simple two sample $t$-test(s). What are the conclussions now?</li>
<li>What would be the corresponding simultaneous confidence intervales for the elements
of the mean vector in this example?</li>
</ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<br><br>


The correponding simultaneous confidence intervals for all possible linear combinations $\boldsymbol{a}^\top \boldsymbol{\mu}$ of the mean vector $\boldsymbol{\mu} \in \mathbb{R}^{p}$ are given as

$$
P\Big(\forall \boldsymbol{a} \in \mathbb{R}^{p};~ \boldsymbol{a}^\top \boldsymbol{\mu} \in \big(
\boldsymbol{a}^\top\overline{\boldsymbol{X}} - \sqrt{K_{\alpha} \boldsymbol{a}^\top \mathcal{S} \boldsymbol{a}}, \boldsymbol{a}^\top\overline{\boldsymbol{X}} + \sqrt{K_{\alpha} \boldsymbol{a}^\top \mathcal{S} \boldsymbol{a}} \big)
\big)\Big) = 1 - \alpha,
$$
where $K_{\alpha}$ is the corresponding quantile obtained from the Fisher $F$ distribution in the form $K_{\alpha} = \frac{p}{n - p} F_{p, n - p}(1 - \alpha)$ and $\mathcal{S}$ is the sample variance-covariance matrix.

```{r}
Fquant <- (9/(20 - 9)) * qf(1 - 0.05, 9, 11)
meanVec <- apply(container.df[,2:10], 2, mean)
LB <- diag(9) %*% meanVec - sqrt(Fquant * diag(diag(9) %*% cov(container.df[,2:10])  %*% diag(9)))
UB <- diag(9) %*% meanVec + sqrt(Fquant * diag(diag(9) %*% cov(container.df[,2:10])  %*% diag(9)))

confInt <- data.frame(LB, UB, meanVec,  row.names = names(container.df)[2:10])
names(confInt) <- c("Lower Boundary", "Upper Boundary", "Mean Estimate")

confInt
```


<br><br>

#### <font color = "#880000"><b>To Do by Yourself</b></font>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<ul>
<li>Consider some other linear combination for the simultaneous confidence intervals and define the necessary linear combinations in order to obtain all pair-wise mean comparissons 
and their corresponding $95~\%$ confidence (simultaneous) intervals.</li><br>

<li>What approach should be used to construct a condidence region (elipsoid)?</li>
</ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<br><br>



#### <font color = "#880000"><b>3. Wilk's Lambda Distribution</b></font>
The Wilks's lambda distribution is defined from two independent Wishart distributed random matrices. The Wilk's lambda distribution is a multivariate 
analogue for the univariate Fisher's F distribution and it is used for inference on two variance-covariace matrices (in a similar way as the Fisher's F distribution 
is used to perform an inference about two variacnce parameters). 

For two independent random matrices 
$$
\mathbb{A} \sim W_{p}(\mathbb{I}, n) \quad \textrm{and} \quad \mathbb{B} \sim W_{p}(\mathbb{I}, m) 
$$

we define a random variable $\frac{|\mathbb{A}|}{|\mathbb{A} + \mathbb{B}|}$ which is said to follow the Wilk's Lambda distribution $\Lambda(p, n, m)$. 
This distribution commnonly appears in the context of likelihood-ratio tests where $n$ is typically the error degrees of freedom, and $m$ is the hypothesis degrees of freedom, so that $n+m$ is the total degrees of freedom.

In a similar way as we perform the analysis of variance for univariate linear regression (ANOVA) we can use a multivariate approach -- so called, multivariate analysis of variance (MANOVA). 
The corresponding inference in MANOVA ($p$-values of the tests) are based on the Wilk's lambda distribution (see the standard R command `manova()` for more details). 



<br><br>

#### <font color = "#880000"><b>To Do by Yourself</b></font>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<ul>
<li>Use some simulations to obtain some usefull insight into the Wilks Lambda distribution.<br> 
(for instance, use the following R code and try to modify it to obtain more complex results)<br>

```{r, fig.width = 10, fig.height=5}
set.seed(1234)
I <- cbind(c(1,0,0), c(0,1,0), c(0,0,1))
A <- rWishart(100, df = 10, I)
B <- rWishart(100, df = 4, I)

WilksLambda <- apply(A, 3, det)/apply(A + B, 3, det)

hist(WilksLambda, col = "lightblue", freq = F, xlab = "Wilk's Lambda", ylab = "Density", ylim = c(0, 4.5))
lines(density(WilksLambda), col = "red", lwd = 2)
```

</li>

<li>Use the univariate case (Wishart's distributions for $p = 1$ and verify (empirically), that 
the Wilk's lambda distribution corresponds with the Fisher's F distribution. </li>
</ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<br><br>

<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: fifth lab session / 25.03.2019)

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 you can assume at least two differen groups. Calculate the estimates for the mean vector in each group 
and apply the statistical test, that both mean vectors are the same. </li>

<li>Choose one group only, and define a reasonable null vector ($\boldsymbol{\mu}_{0}$) for the true mean vector ($\boldsymbol{\mu}$) 
and test the null hypothesis that $\boldsymbol{\mu}$ is equal to $\boldsymbol{\mu}_{0}$. Provide the corresponding $95~\%$ confidence region and interpret the results. </li>
</ul>




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