---
title: "P3.2 - Ejercicio 10"
author: "Paul Rognon"
date: "06/10/2019"
output:
  html_document: default
editor_options:
  chunk_output_type: console
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE,warning = FALSE)
remove(list=ls())
```


### Hipótesis a estudiar
**H01**: Los dos perfiles son similares (paralelos)

$$\forall k=1...p-1,\hspace{0.3cm}\mu_{tnt,k+1}-\mu_{tnt,k}=\mu_{c,k+1}-\mu_{c,k}$$

**H02**: Los dos perfiles están al mismo nivel
$$\frac{1}{p}\sum_{j=1}^{p}{\mu_{tnt,j}}=\frac{1}{p}\sum_{j=1}^{p}{\mu_{c,j}}$$

**H03**: No hay diferencia entre las medias de los tests (no hay efecto principal por test)
$$\frac{\mu_{tnt,1}+\mu_{c,1}}{2}=\frac{\mu_{tnt,2}+\mu_{c,2}}{2}=...=\frac{\mu_{tnt,p}+\mu_{c,p}}{2}$$

El gráfico de los perfiles es:
```{r,include=FALSE}
remove(list=ls())
sample1<-read.table("ej4pra3-c.txt",row.names = 1)
sample2<-read.table("ej4pra3-tnt.txt",row.names = 1)

#b# Grafico de perfiles
sample11<-sample1
sample22<-sample2
sample11$grupo<-rep("C",nrow(sample11))
sample22$grupo<-rep("TNT",nrow(sample22))
data2<-rbind(sample11,sample22)

library(reshape2)
library(dplyr)
library(ggplot2)

perfiles<-aggregate(data2[,1:ncol(sample1)], list(data2$grupo), mean)
data3 <- melt(perfiles)
plo<-ggplot(data3) +
  geom_line(aes(variable, value, group = Group.1,color = Group.1))

#c# Test H01
n1<-nrow(sample1)
n2<-nrow(sample2)
x_barra1<-colMeans(sample1)
x_barra2<-colMeans(sample2)
p<-qr(sample1)$rank

Q_1<-matrix(rep(0,p*p),ncol=p)
for (j  in 1:nrow(sample1)){
  Q_1<-Q_1+(t(sample1[j,])-x_barra1)%*%t(t(sample1[j,])-x_barra1)
}
S_1<-Q_1/(n1-1)

Q_2<-matrix(rep(0,p*p),ncol=p)
for (j  in 1:nrow(sample2)){
  Q_2<-Q_2+(t(sample2[j,])-x_barra2)%*%t(t(sample2[j,])-x_barra2)
}
S_2<-Q_2/(n2-1)

S<-((n1-1)*S_1+(n2-1)*S_2)/(n1+n2-2)

C<-rbind(c(1,-1,0,0,0,0,0,0,0,0),
         c(0,1,-1,0,0,0,0,0,0,0),
         c(0,0,1,-1,0,0,0,0,0,0),
         c(0,0,0,1,-1,0,0,0,0,0),
         c(0,0,0,0,1,-1,0,0,0,0),
         c(0,0,0,0,0,1,-1,0,0,0),
         c(0,0,0,0,0,0,1,-1,0,0),
         c(0,0,0,0,0,0,0,1,-1,0),
         c(0,0,0,0,0,0,0,0,1,-1)
)

q<-qr(C)$rank
mu_0<-rep(0,q)


T0_2_y<-t((C%*%(x_barra1-x_barra2))-mu_0)%*%solve(C%*%S%*%t(C))%*%(C%*%((x_barra1-x_barra2))-mu_0)*n1*n2/(n1+n2)
F0<-(n1+n2-2-q+1)*T0_2_y/(q*(n1+n2-2))
f_alpha<-qf(0.05,q,n1+n2-2-q+1,lower.tail = FALSE)
p_valor_H01<-pf(F0,q,n1+n2-2-q+1,lower.tail = FALSE)
#print(p_valor)
#p-valor de 0.5537363, no rechazo H01

#c# Tests de H03
#Testear H03 es como testar que la suma ponderada de tnt y c tiene como vector de media un multiple del vector de unos.
n1<-nrow(sample1)
n2<-nrow(sample2)

x_barra1<-colMeans(sample1)
x_barra2<-colMeans(sample2)
p<-qr(sample1)$rank


z_barra<-(n1*x_barra1+n2*x_barra2)/(n1+n2)

Q_1<-matrix(rep(0,p*p),ncol=p)
for (j  in 1:nrow(sample1)){
  Q_1<-Q_1+(t(sample1[j,])-x_barra1)%*%t(t(sample1[j,])-x_barra1)
}
S_1<-Q_1/(n1-1)

Q_2<-matrix(rep(0,p*p),ncol=p)
for (j  in 1:nrow(sample2)){
  Q_2<-Q_2+(t(sample2[j,])-x_barra2)%*%t(t(sample2[j,])-x_barra2)
}
S_2<-Q_2/(n2-1)

S<-((n1-1)*S_1+(n2-1)*S_2)/(n1+n2-2)

C<-rbind(c(1,-1,0,0,0,0,0,0,0,0),
         c(0,1,-1,0,0,0,0,0,0,0),
         c(0,0,1,-1,0,0,0,0,0,0),
         c(0,0,0,1,-1,0,0,0,0,0),
         c(0,0,0,0,1,-1,0,0,0,0),
         c(0,0,0,0,0,1,-1,0,0,0),
         c(0,0,0,0,0,0,1,-1,0,0),
         c(0,0,0,0,0,0,0,1,-1,0),
         c(0,0,0,0,0,0,0,0,1,-1)
         )

q<-qr(C)$rank
mu_0<-rep(0,q)


T0_2_y<-t((C%*%z_barra)-mu_0)%*%solve(C%*%S%*%t(C))%*%(C%*%z_barra-mu_0)*(n1+n2)
F0<-(n1+n2-2-q+1)*T0_2_y/(q*(n1+n2-2))
f_alpha<-qf(0.05,q,n1+n2-2-q+1,lower.tail = FALSE)
p_valor_H03<-pf(F0,q,n1+n2-2-q+1,lower.tail = FALSE)
#print(p_valor)
#p-valor de 1.943049e-09, Rechazo H03
```
```{r,echo=FALSE}
plo
```


### Resumen del Ejercicio P3.1 - 4
No se había rechazado H01, con un p-valor de `r round(p_valor_H01,4)`.
Entonces se había estudiado la hipótesis H03 que se había rechazado con un p-valor de $\approx$ `r round(p_valor_H03,4)`


### Ejercicio P3.2 - 10 Punto a

#### *¿es necesario suponer que las matrices de covarianza en las dos poblaciones sean las mismas o se puede hacer un supuesto más débil?*

Para el test de igualdad de medias entre dos poblaciones, el supuesto de igualdad de las matrices de covarianza se puede relajar sí:

1. si las dos poblaciones son del mismo tamaño: 
    + $n_1=n_2=n$ y $n$ suficientemente grande, entonces el estadístico $T_0^2$ es asymptoticamente insensible a $\Sigma_1 \neq \Sigma_2$.
    + se puede correr un test sobre la media de la variable $y=x_1-x_2$ que tiene una distribución $N_p(\mu_1-\mu_2,\Sigma_1 + \Sigma_2)$. Se testea $H_0 : \mu_y=0$. El estadístico resultante tiene distribución $T_{p,n-1}^2$ mientras que en la solución previa, el estadístico se basa en $2(n-1)$ grados de libertad. 
2. si $n_1\approx n_2$ y son muy grandes, $\Sigma_1 \neq \Sigma_2$ no impacta el nivel de significación y el poder del test.

Sino hay que usar una solución al problema de Fisher-Behrens, con la de James, Yao o un bootstrap paramétrico.

En nuestro caso tenemos $n_1=n_2=n=27$. Si consideramos que 27 es suficientemente grande y buscamos aplicar la solución 1a. 
Primero notamos que en el caso  $n_1=n_2=n$ y $\Sigma_1=\Sigma_2$, bajo $H_0$ la distribución del estadístico de Hotelling para dos muestras tiende a una $\chi^2_p$:

$T_0^2= \frac{(n-2)p}{n-p-1}F_{p,n-p-1}$

cuando $n\longrightarrow+\infty$

  $$\frac{(n-2)p}{n-p-1}\longrightarrow p$$

  $$T_0^2= pF_{p,n-p-1 \approx }\chi^2_p$$

Cuando $n_1=n_2=n$ es suficientemente grande se puede suponer $S_1=\Sigma_1$ y $S_2=\Sigma_2$. Entonces bajo, $H_0$ el estadístico de de Hotelling para dos muestras:

$$T_0^2=\frac{n_1n_2}{n_1+n_2}(\bar{x_1}-\bar{x_2})^T(\frac{(n_1-1)S_1+(n_2-2)S_2}{n_1+n_2-2})^{-1}(\bar{x_1}-\bar{x_2})$$

$$T_0^2=\frac{n_1+n_2-2}{n_1+n_2}(\bar{x_1}-\bar{x_2})^T(\frac{n_1-1}{n_1n_2}S_1+\frac{n_2-1}{n_1n_2}{S_2})^{-1}(\bar{x_1}-\bar{x_2})$$

$$T_0^2=(\bar{x_1}-\bar{x_2})^{T}(\frac{1}{n_2} \Sigma_1+\frac{1}{n_1} \Sigma_2)^{-1}(\bar{x_1}-\bar{x_2})$$

Pero:

$\bar{x_1}-\bar{x_2} \sim N_p(0,\frac{1}{n_1}\Sigma_1+\frac{1}{n_2}\Sigma_2)$

Entonces:

$(\bar{x_1}-\bar{x_2})^T(\frac{1}{n_1}\Sigma_1+\frac{1}{n_2}\Sigma_2)^{-1}(\bar{x_1}-\bar{x_2}) \sim \chi^2_p$

Cuando $n_1=n_2=n$ suficientemente grande, el estadístico T tiene la distribución asintótica correcta.

#### *Exprese matemáticamente la hipótesis sobre la igualdad de covarianzas que sea necesaria en cada caso*

En los tests para **H01**, **H02** y **H03**, se testea la igualdad de medias de transformaciones de $x_1$ y $x_2$.

**Para H01**:  
Definimos $C_1$:

$$
C_1=\left(\begin{array}{cccccccccc} 
1&-1&0&0&0&0&0&0&0&0\\
0&1&-1&0&0&0&0&0&0&0\\
0&0&1&-1&0&0&0&0&0&0\\
0&0&0&1&-1&0&0&0&0&0\\
0&0&0&0&1&-1&0&0&0&0\\
0&0&0&0&0&1&-1&0&0&0\\
0&0&0&0&0&0&1&-1&0&0\\
0&0&0&0&0&0&0&1&-1&0\\
0&0&0&0&0&0&0&0&1&-1
\end{array}\right)
$$

$y_1=C_1x_1, \hspace{0.5cm}y_1\sim N_p(C_1\mu_1,C_1\Sigma_1C_1^T)$

$y_2=C_1x_2, \hspace{0.5cm}y_2\sim N_p(C_1\mu_2,C_1\Sigma_2C_1^T)$

**H01** equivale es testear: $\mu_{y1}=\mu_{y2}$

La hipótesis sobre covarianzas necesaria es entonces
$$\Sigma_{y1}=\Sigma_{y2}\iff C_1\Sigma_1C_1^T=C_1\Sigma_2C_1^T$$

**Para H02**:  
Definimos $C_2$:

$$
\begin{align}
    C_2 &= (1,1,1,1,1,1,1,1,1,1)
\end{align}
$$

$y_1=C_2x_1, \hspace{0.5cm}y_1\sim N_p(C_2\mu_1,C_2\Sigma_1C_2^T)$

$y_2=C_2x_2, \hspace{0.5cm}y_2\sim N_p(C_2\mu_2,C_2\Sigma_2C_2^T)$

**H02** equivale es testear: $\mu_{y1}=\mu_{y2}$

La hipótesis sobre covarianzas necesaria es entonces $$\Sigma_{y1}=\Sigma_{y2}\iff C_2\Sigma_1C_2^T=C_2\Sigma_2C_2^T\iff\sum_{i,j=1}^{p}{\sigma_{1,i,j}}=\sum_{i,j=1}^{p}{\sigma_{2,i,j}}$$


**Para H03**:  
Definimos $C_1$:
$$
C_1=\left(\begin{array}{cccccccccc} 
1&-1&0&0&0&0&0&0&0&0\\
0&1&-1&0&0&0&0&0&0&0\\
0&0&1&-1&0&0&0&0&0&0\\
0&0&0&1&-1&0&0&0&0&0\\
0&0&0&0&1&-1&0&0&0&0\\
0&0&0&0&0&1&-1&0&0&0\\
0&0&0&0&0&0&1&-1&0&0\\
0&0&0&0&0&0&0&1&-1&0\\
0&0&0&0&0&0&0&0&1&-1
\end{array}\right)
$$

$y_1=C_1x_1, \hspace{0.5cm}y_1\sim N_p(C_1\mu_1,C_1\Sigma_1C_1^T)$

$y_2=-C_1x_2, \hspace{0.5cm}y_2\sim N_p(-C_1\mu_2,C_1\Sigma_2C_1^T)$

**H03** equivale es testear: $\mu_{y1}=\mu_{y2}$

La hipótesis sobre covarianzas necesaria es entonces 
$$\Sigma_{y1}=\Sigma_{y2}\iff C_1\Sigma_1C_1^T=C_1\Sigma_2C_1^T$$

### Ejercicio P3.2 - 10 Punto b

#### *Para cada una de las situaciones de interés estudie si el supuesto de igualdad de covarianzas es razonable.*

**Para H01 y H03**  
Para facilitar el cáclulo en R usamos que como $n_1=n_2=n$:

$$\gamma^*=\frac{|\frac{Q_1}{n_1}|^{\frac{n_1}{2}}|\frac{Q_2}{n_2}|^{\frac{n_2}{2}}}{|\frac{U}{n_1+n_2}|^{\frac{n_1+n_2}{2}}}=\frac{|\frac{Q_1}{n}|^{\frac{n}{2}}|\frac{Q_2}{n}|^{\frac{n}{2}}}{|\frac{U}{2n}|^{\frac{2n}{2}}}=\left(\frac{|\frac{Q_1}{n}||\frac{Q_2}{n}|}{|\frac{U}{2n}|}\right)^{\frac{n}{2}}\frac{1}{|\frac{U}{2n}|^{\frac{n}{2}}}$$
```{r}
remove(list=ls())
sample1<-read.table("ej4pra3-c.txt",row.names = 1)
sample2<-read.table("ej4pra3-tnt.txt",row.names = 1)

C<-rbind(c(1,-1,0,0,0,0,0,0,0,0),
         c(0,1,-1,0,0,0,0,0,0,0),
         c(0,0,1,-1,0,0,0,0,0,0),
         c(0,0,0,1,-1,0,0,0,0,0),
         c(0,0,0,0,1,-1,0,0,0,0),
         c(0,0,0,0,0,1,-1,0,0,0),
         c(0,0,0,0,0,0,1,-1,0,0),
         c(0,0,0,0,0,0,0,1,-1,0),
         c(0,0,0,0,0,0,0,0,1,-1)
)

sample1<-t(C%*%t(sample1))
sample2<-t(C%*%t(sample2))
# Test hipotesisis de covarianzas para H01
n1<-nrow(sample1)
n2<-nrow(sample2)
x_barra1<-colMeans(sample1)
x_barra2<-colMeans(sample2)
p<-qr(sample1)$rank

Q_1<-matrix(rep(0,p*p),ncol=p)
for (j  in 1:nrow(sample1)){
  Q_1<-Q_1+(sample1[j,]-x_barra1)%*%t(sample1[j,]-x_barra1)
}


Q_2<-matrix(rep(0,p*p),ncol=p)
for (j  in 1:nrow(sample2)){
  Q_2<-Q_2+(sample2[j,]-x_barra2)%*%t(sample2[j,]-x_barra2)
}


U<-Q_1+Q_2

k<-2

n<-n1+n2
estadis<-log(((det(Q_1/n1)*det(Q_2/n2)/det(U/(2*n1)))^(n1/2))/(det(U/(2*n1))^(n1/2)))*(-2)

nu<-k*p*(p+1)/2
nu_1<-p*(p+1)/2

chi_alpha<-qchisq(0.01,nu-nu_1,lower.tail = FALSE)
p_valor_H01<-pchisq(estadis,nu-nu_1,lower.tail = FALSE)
```
Obtenemos un p-valor de `r round(p_valor_H01,4)`, entonces no rechazamos al 1%.

**Para H02**
```{r}
remove(list=ls())
sample1<-read.table("ej4pra3-c.txt",row.names = 1)
sample2<-read.table("ej4pra3-tnt.txt",row.names = 1)

C<-rep(1,10)
sample1<-t(C%*%t(sample1))
sample2<-t(C%*%t(sample2))
# Test hipotesisis de covarianzas para H01
n1<-nrow(sample1)
n2<-nrow(sample2)
x_barra1<-colMeans(sample1)
x_barra2<-colMeans(sample2)
p<-qr(sample1)$rank

Q_1<-matrix(rep(0,p*p),ncol=p)
for (j  in 1:nrow(sample1)){
  Q_1<-Q_1+(sample1[j,]-x_barra1)%*%t(sample1[j,]-x_barra1)
}


Q_2<-matrix(rep(0,p*p),ncol=p)
for (j  in 1:nrow(sample2)){
  Q_2<-Q_2+(sample2[j,]-x_barra2)%*%t(sample2[j,]-x_barra2)
}


U<-Q_1+Q_2

k<-2

n<-n1+n2
estadis<-log(((det(Q_1/n1)*det(Q_2/n2)/det(U/(2*n1)))^(n1/2))/(det(U/(2*n1))^(n1/2)))*(-2)

nu<-k*p*(p+1)/2
nu_1<-p*(p+1)/2

chi_alpha<-qchisq(0.01,nu-nu_1,lower.tail = FALSE)
p_valor_H02<-pchisq(estadis,nu-nu_1,lower.tail = FALSE)
```
Obtenemos un p-valor de `r round(p_valor_H02,4)`, entonces no rechazamos al 1%.

#### Verificamos el supuesto de normalidad
```{r, include=FALSE}
remove(list=ls())
sample1<-read.table("ej4pra3-c.txt",row.names = 1)
sample2<-read.table("ej4pra3-tnt.txt",row.names = 1)
library(mvShapiroTest)
s_res1<-mvShapiro.Test(as.matrix(sample1)) 
s_res2<-mvShapiro.Test(as.matrix(sample2))
```
Corriendo un test de Shapiro multivariado sobre la muestra bruta, obtenemos un p-valor de `r round(s_res1$p.value,4)` para la muestra C. Obtenemos un p-valor de `r round(s_res2$p.value,4)` para la muestra TNT. Rechazo la normalidad al 5% para la muestra C y rechazo al 20% para la muestra TNT.

En realidad lo que importa es la normalidad de las variables transformadas que son las que testeamos efectivamente. Pasamos entonces a testear normalidad para las variables transformadas.

**Para H01 y H03**  
```{r}
remove(list=ls())
sample1<-read.table("ej4pra3-c.txt",row.names = 1)
sample2<-read.table("ej4pra3-tnt.txt",row.names = 1)
C<-rbind(c(1,-1,0,0,0,0,0,0,0,0),
         c(0,1,-1,0,0,0,0,0,0,0),
         c(0,0,1,-1,0,0,0,0,0,0),
         c(0,0,0,1,-1,0,0,0,0,0),
         c(0,0,0,0,1,-1,0,0,0,0),
         c(0,0,0,0,0,1,-1,0,0,0),
         c(0,0,0,0,0,0,1,-1,0,0),
         c(0,0,0,0,0,0,0,1,-1,0),
         c(0,0,0,0,0,0,0,0,1,-1)
)

sample1<-t(C%*%t(sample1))
sample2<-t(C%*%t(sample2))
library(mvShapiroTest)
s_res1<-mvShapiro.Test(as.matrix(sample1)) 
s_res2<-mvShapiro.Test(as.matrix(sample2))
```
Corriendo un test de Shapiro multivariado, obtenemos un p-valor de `r round(s_res1$p.value,4)` para la muestra C. Obtenemos un p-valor de `r round(s_res2$p.value,4)` para la muestra TNT. En ambos casos no rechazo la normalidad al 20%.

```{r}
remove(list=ls())
sample1<-read.table("ej4pra3-c.txt",row.names = 1)
sample2<-read.table("ej4pra3-tnt.txt",row.names = 1)

C<-rep(1,10)
sample1<-(t(C%*%t(sample1)))[,1]
sample2<-(t(C%*%t(sample2)))[,1]

#Puedo correr un test univariado porque son vectores
s_res1<-shapiro.test(sample1)
s_res2<-shapiro.test(sample2)
```
Corriendo un test de Shapiro multivariado, obtenemos un p-valor de `r round(s_res1$p.value,4)` para la muestra C. Obtenemos un p-valor de `r round(s_res2$p.value,4)` para la muestra TNT. En ambos casos no rechazo la normalidad al 20%.

Era válido aplicar el test de igualdad de matrices de covarianzas para H01, H02 y H03.


### TODO LOS QUE ESTÁ ABAJO NO SIRVE
No sirve porque no rechazamos normalidad para las variables que testeamos efectivamente. Lo dejo por si el código les sirve.

Calculamos un estimador robusto de la matriz de covarianza de la muestra C usando covRob de RobStatTM. Realizamos un QQ-plot de las distancias de Mahalanobis robustas obtenidas contra una distribución Chi-cuadrado con 10 grados de libertad.
```{r, include=FALSE}
remove(list=ls())
sample1<-read.table("ej4pra3-c.txt",row.names = 1)
sample2<-read.table("ej4pra3-tnt.txt",row.names = 1)
library(RobStatTM)
library(car)
(res<-initPP(sample1))
sum(res$idx)#initPP da 7 outliers

cov_rob_res<-covRob(sample1)
p<-qr(sample1)$rank
chisq_p_sim<-rchisq(150,p)
windows()
par(mfrow=c(1,2),oma=c(0,0,2,0))
plot_qq<-qqplot(chisq_p_sim,cov_rob_res$dist,xlim=c(0,50),ylim=c(0,50),xlab="Chi-cuadrado",ylab="Distancias")
Boxplot(cov_rob_res$dist)
title("Muestra C",outer=T)
#identify(plot_qq$x,plot_qq$y,n=2)
savePlot("outliers in sample1","png")
dev.off()


(res<-initPP(sample2))
sum(res$idx)#initPP da 6 outliers

p<-qr(sample2)$rank
chisq_p_sim<-rchisq(150,p)
windows()
par(mfrow=c(1,2),oma=c(0,0,2,0))
plot_qq<-qqplot(chisq_p_sim,res$disma)
Boxplot(res$disma)
title("Muestra TNT",outer=T)
#identify(plot_qq$x,plot_qq$y,n=6)
savePlot("outliers in sample2","png")
dev.off()
```

```{r,echo=FALSE}
knitr::include_graphics("outliers in sample1.png")
knitr::include_graphics("outliers in sample2.png")
#![ ]("\outliers in sample1.png")
```
Quitamos los outliers identificados. Volvemos a graficar QQ-Plot y Boxplot hasta obtener los gráficos siguientes donde quedan un sólo outlier en cada muestra
```{r,include=FALSE}
remove(list=ls())
sample1<-read.table("ej4pra3-c.txt",row.names = 1)
sample1<-sample1[-c(18,19,26),]
sample1<-sample1[-c(10,11,7,13),]
sample2<-read.table("ej4pra3-tnt.txt",row.names = 1)
sample2<-sample2[-c(5,6,15,20,22,23),]
sample2<-sample2[-c(18),]

(res<-initPP(sample1))
cov_rob_res<-covRob(sample1)
p<-qr(sample1)$rank
chisq_p_sim<-rchisq(150,p)
windows()
par(mfrow=c(1,2),oma=c(0,0,2,0))
plot_qq<-qqplot(chisq_p_sim,res$disma,xlab="Chi-cuadrado",ylab="Distancias")
i<-which.max(res$disma)
text(plot_qq$x[length(plot_qq$y)],plot_qq$y[length(plot_qq$y)]-0.5,as.character(i))
Boxplot(res$disma)
title("Muestra C",outer=T)
#identify(plot_qq$x,plot_qq$y,n=3)
savePlot("outliers in sample11","png")
dev.off()
sample1<-sample1[-c(2),]


(res<-initPP(sample2))
sum(res$idx)#initPP da 6 outliers
cov_rob_res<-covRob(sample2)

p<-qr(sample2)$rank
chisq_p_sim<-rchisq(150,p)
windows()
par(mfrow=c(1,2),oma=c(0,0,2,0))
plot_qq<-qqplot(chisq_p_sim,cov_rob_res$dist)
i<-which.max(cov_rob_res$dist)
text(plot_qq$x[length(plot_qq$y)],plot_qq$y[length(plot_qq$y)]-0.5,as.character(i))
Boxplot(cov_rob_res$dist)
title("Muestra TNT",outer=T)
#identify(plot_qq$x,plot_qq$y,n=2)
savePlot("outliers in sample22","png")
dev.off()
sample2<-sample2[-c(9),]
```

```{r,echo=FALSE}
knitr::include_graphics("outliers in sample11.png")
knitr::include_graphics("outliers in sample22.png")
#![ ]("\outliers in sample1.png")

s_res1<-mvShapiro.Test(as.matrix(sample1))
s_res2<-mvShapiro.Test(as.matrix(sample2))
```
Corriendo un test de Shapiro multivariado de nuevo, obtenemos un p-valor de `r round(s_res1$p.value,4)` para la muestra C. Obtenemos un p-valor de `r round(s_res2$p.value,4)` para la muestra TNT. No rechazamos normalidad al 20%

#### Vuelvo a correr los tests sobre covarianzas sin los outliers.

**Para H01 y H03**
```{r}
remove(list=ls())
sample1<-read.table("ej4pra3-c.txt",row.names = 1)
sample1<-sample1[-c(18,19,26),]
sample1<-sample1[-c(10,11,7,13),]
sample1<-sample1[-c(2),]
sample2<-read.table("ej4pra3-tnt.txt",row.names = 1)
sample2<-sample2[-c(5,6,15,20,22,23),]
sample2<-sample2[-c(18),]
sample2<-sample2[-c(9),]

C<-rbind(c(1,-1,0,0,0,0,0,0,0,0),
         c(0,1,-1,0,0,0,0,0,0,0),
         c(0,0,1,-1,0,0,0,0,0,0),
         c(0,0,0,1,-1,0,0,0,0,0),
         c(0,0,0,0,1,-1,0,0,0,0),
         c(0,0,0,0,0,1,-1,0,0,0),
         c(0,0,0,0,0,0,1,-1,0,0),
         c(0,0,0,0,0,0,0,1,-1,0),
         c(0,0,0,0,0,0,0,0,1,-1)
)

sample1<-t(C%*%t(sample1))
sample2<-t(C%*%t(sample2))
# Test hipotesisis de covarianzas para H01
n1<-nrow(sample1)
n2<-nrow(sample2)
x_barra1<-colMeans(sample1)
x_barra2<-colMeans(sample2)
p<-qr(sample1)$rank

Q_1<-matrix(rep(0,p*p),ncol=p)
for (j  in 1:nrow(sample1)){
  Q_1<-Q_1+(sample1[j,]-x_barra1)%*%t(sample1[j,]-x_barra1)
}


Q_2<-matrix(rep(0,p*p),ncol=p)
for (j  in 1:nrow(sample2)){
  Q_2<-Q_2+(sample2[j,]-x_barra2)%*%t(sample2[j,]-x_barra2)
}


U<-Q_1+Q_2

k<-2

nm<-min(n1,n2)
n<-n1+n2
estadis<-log(((det(Q_1/n1)*det(Q_2/n2)/det(U/n))^(nm/2))*((det(Q_1/n1)^((n1-nm)/2))*(det(Q_2/n2)^((n2-nm)/2))/(det(U/n)^((n-nm)/2))))*(-2)

nu<-k*p*(p+1)/2
nu_1<-p*(p+1)/2

chi_alpha<-qchisq(0.01,nu-nu_1,lower.tail = FALSE)
p_valor_H01<-pchisq(estadis,nu-nu_1,lower.tail = FALSE)
p_valor_H01
```


Obtenemos un p-valor de approx `r round(p_valor_H01,4)`, entonces rechazamos al 1%.

**Para H02**
```{r}
remove(list=ls())
sample1<-read.table("ej4pra3-c.txt",row.names = 1)
sample1<-sample1[-c(18,19,26),]
sample1<-sample1[-c(10,11,7,13),]
sample1<-sample1[-c(2),]
sample2<-read.table("ej4pra3-tnt.txt",row.names = 1)
sample2<-sample2[-c(5,6,15,20,22,23),]
sample2<-sample2[-c(18),]
sample2<-sample2[-c(9),]

C<-rep(1,10)
sample1<-t(C%*%t(sample1))
sample2<-t(C%*%t(sample2))
# Test hipotesisis de covarianzas para H01
n1<-nrow(sample1)
n2<-nrow(sample2)
x_barra1<-colMeans(sample1)
x_barra2<-colMeans(sample2)
p<-qr(sample1)$rank

Q_1<-matrix(rep(0,p*p),ncol=p)
for (j  in 1:nrow(sample1)){
  Q_1<-Q_1+(sample1[j,]-x_barra1)%*%t(sample1[j,]-x_barra1)
}


Q_2<-matrix(rep(0,p*p),ncol=p)
for (j  in 1:nrow(sample2)){
  Q_2<-Q_2+(sample2[j,]-x_barra2)%*%t(sample2[j,]-x_barra2)
}


U<-Q_1+Q_2

k<-2

nm<-min(n1,n2)
n<-n1+n2
estadis<-log(((det(Q_1/n1)*det(Q_2/n2)/det(U/n))^(nm/2))*((det(Q_1/n1)^((n1-nm)/2))*(det(Q_2/n2)^((n2-nm)/2))/(det(U/n)^((n-nm)/2))))*(-2)

nu<-k*p*(p+1)/2
nu_1<-p*(p+1)/2

chi_alpha<-qchisq(0.01,nu-nu_1,lower.tail = FALSE)

p_valor_H02<-pchisq(estadis,nu-nu_1,lower.tail = FALSE)
p_valor_H02
```
Obtenemos un p-valor de approx.`r round(p_valor_H02,4)`, entonces no rechazamos al 1%.
