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