---
title: "Practica 1 - Analisis Multivariado I"
author: "Ciro Carvallo"
date: "Septiembre 5, 2019"
output:
  html_document: default
  word_document: default
  pdf_document: default
---

```{r}

densidad_hotelling <- function(y,p,m,ncp = FALSE) {
  if (ncp == FALSE){
    return ((m- p +1)*df(y*((m- p +1)/(p*m)), p, m - p + 1, log = FALSE)/(p*m))
  }
  else{
    return ((m- p +1)*df(y*((m- p +1)/(p*m)), p, m - p + 1, ncp, log = FALSE)/(p*m))
  }
}

```


# Practica 1 - Parte 2
## Ejercicio 1.2 2 

```{r}
  library(MASS)
  library(rWishart)
  library(ggplot2)
```

### Parte A
```{r}

  set.seed(42)
  NITER = 1000
  vector_hotelling = rep(NA,NITER)
  for (i in 1:NITER){
      mu = c(0,0,0)
      Sigma = cbind(c(5,0,2),c(0,5,1),c(2,1,2))
      mvnormal = mvrnorm(n = 1, mu, Sigma) # n = sample size
      wishart = rWishart(n = 1,df = 20,Sigma)
      hotelling = 20*t(mvnormal)%*%solve(wishart[,,1])%*%mvnormal
      vector_hotelling[i] = hotelling

  }
  #hist(vector_hotelling)
  #print(density(vector_hotelling))
  df <-data.frame(hotelling = vector_hotelling)
  ggplot(df, aes(x=hotelling)) + geom_histogram(aes(y = ..density..),
  colour = "blue", fill = "white") +
  geom_density(alpha = .2, fill="#FF6655")
  y =seq(from = 0, to = 30, by = 0.05)
  densidad = densidad_hotelling(y,3,20)
  df2 = data.frame(y,densidad )
  ggplot(df, aes(x=hotelling)) + geom_histogram(aes(y = ..density..),
  colour = "blue", fill = "white") +
  geom_line(data = df2,aes(x = y , y = densidad)) 
  #print(df(y*(20-3+1)/(20*3),3,20-3+1)*(20-3+1)/(20*3))
  
```


## Parte B

```{r}
  set.seed(42)
  n = 3
  NITER = 1000
  vector_hotelling = rep(NA,NITER)
  for (i in 1:NITER){
      mu = c(-1,2,-3)
      Sigma = cbind(c(5,0,2),c(0,5,1),c(2,1,2))
      mvnormal = mvrnorm(n = 1, mu, Sigma) # n = sample size
      wishart = rWishart(n = 1,df = 20,Sigma)
      hotelling = 20*t(mvnormal)%*%solve(wishart[,,1])%*%mvnormal
      vector_hotelling[i] = hotelling
  }
  #hist(vector_hotelling)
  #print(density(vector_hotelling))
  df <-data.frame(hotelling = vector_hotelling)
  ggplot(df, aes(x=hotelling)) + geom_histogram(aes(y = ..density..),
  colour = "blue", fill = "white") +
  geom_density(alpha = .2, fill="#FF6655")
  
  ncp = t(mu)%*%solve(Sigma)%*%mu
  
  y =seq(from = 0, to = 80, by = 0.05)
  densidad = densidad_hotelling(y,3,20,ncp)
  df2 = data.frame(y,densidad )
  ggplot(df, aes(x=hotelling)) + geom_histogram(aes(y = ..density..),
  colour = "blue", fill = "white") +
  geom_line(data = df2,aes(x = y , y = densidad))

```

