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

  library(MASS)
  library(rWishart)
## Warning: package 'rWishart' was built under R version 3.6.1
## 
## Attaching package: 'rWishart'
## The following object is masked from 'package:stats':
## 
##     rWishart
  library(ggplot2)
## Registered S3 methods overwritten by 'ggplot2':
##   method         from 
##   [.quosures     rlang
##   c.quosures     rlang
##   print.quosures rlang

Parte A

  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")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

  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)) 
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

  #print(df(y*(20-3+1)/(20*3),3,20-3+1)*(20-3+1)/(20*3))

Parte B

  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")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

  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))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.