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