data(mtcars)
X <- mtcars[,c("mpg","disp","hp","wt")]
#Hipotesis
Ho<-c(20,200,150,3) #rata-rata populasi berdasarkan hipotesis awal
X
##                      mpg  disp  hp    wt
## Mazda RX4           21.0 160.0 110 2.620
## Mazda RX4 Wag       21.0 160.0 110 2.875
## Datsun 710          22.8 108.0  93 2.320
## Hornet 4 Drive      21.4 258.0 110 3.215
## Hornet Sportabout   18.7 360.0 175 3.440
## Valiant             18.1 225.0 105 3.460
## Duster 360          14.3 360.0 245 3.570
## Merc 240D           24.4 146.7  62 3.190
## Merc 230            22.8 140.8  95 3.150
## Merc 280            19.2 167.6 123 3.440
## Merc 280C           17.8 167.6 123 3.440
## Merc 450SE          16.4 275.8 180 4.070
## Merc 450SL          17.3 275.8 180 3.730
## Merc 450SLC         15.2 275.8 180 3.780
## Cadillac Fleetwood  10.4 472.0 205 5.250
## Lincoln Continental 10.4 460.0 215 5.424
## Chrysler Imperial   14.7 440.0 230 5.345
## Fiat 128            32.4  78.7  66 2.200
## Honda Civic         30.4  75.7  52 1.615
## Toyota Corolla      33.9  71.1  65 1.835
## Toyota Corona       21.5 120.1  97 2.465
## Dodge Challenger    15.5 318.0 150 3.520
## AMC Javelin         15.2 304.0 150 3.435
## Camaro Z28          13.3 350.0 245 3.840
## Pontiac Firebird    19.2 400.0 175 3.845
## Fiat X1-9           27.3  79.0  66 1.935
## Porsche 914-2       26.0 120.3  91 2.140
## Lotus Europa        30.4  95.1 113 1.513
## Ford Pantera L      15.8 351.0 264 3.170
## Ferrari Dino        19.7 145.0 175 2.770
## Maserati Bora       15.0 301.0 335 3.570
## Volvo 142E          21.4 121.0 109 2.780
Hotelling <- function(X, mu=0, asymp=FALSE,alfa=0.05){
X <- as.matrix(X)
n <- nrow(X)
p <- ncol(X)
df2 <- n - p
if(df2 < 1L) stop("Need nrow(X) > ncol(X).")
if(length(mu) != p) mu <- rep(mu[1], p)
xbar <- colMeans(X)
S <- cov(X)
T2 <- n * t(xbar - mu) %*% solve(S) %*% (xbar - mu)
#Statistik Uji
Fstat <- T2 / (p * (n-1) / df2)
if(asymp){
pval <- 1 - pchisq(T2, df=p)
} else {
pval <- 1 - pf(Fstat, df1=p, df2=df2)
}
data.frame(T2=as.numeric(T2), Fstat=as.numeric(Fstat),df1=p, df2=df2, p.value=as.numeric(pval),
asymp=asymp, row.names="")
}
Hotelling(X,mu=Ho,asymp =TRUE)
##        T2    Fstat df1 df2   p.value asymp
##  10.78587 2.435519   4  28 0.0290789  TRUE
library("DescTools")
## Warning: package 'DescTools' was built under R version 3.6.3
#test f and chi
HotellingsT2Test(X, mu=Ho, test = 'chi')
## 
##  Hotelling's one sample T2-test
## 
## data:  X
## T.2 = 10.786, df = 4, p-value = 0.02908
## alternative hypothesis: true location is not equal to c(20,200,150,3)
#see dataset X on 2 variable
n <- nrow(X)
p <- ncol(X)
xbar <- colMeans(X)
S <- cov(X)
library(car)
## Loading required package: carData
## Warning: package 'carData' was built under R version 3.6.3
## 
## Attaching package: 'car'
## The following object is masked from 'package:DescTools':
## 
##     Recode
alpha<-0.05
tconst <- sqrt((p/n)*((n-1)/(n-p)) * qf(1-alpha,p,n-p))
id <- c(1,2)
plot(ellipse(center=xbar[id], shape=S[id,id],
radius=tconst, draw=F),type="n", xlab="mpg", ylab="disp")
lines(ellipse(center=xbar[id], shape=S[id,id],
radius=tconst, lwd=3), xlab="mpg", ylab="disp")
points(20,200,col="red", pch=16)
text(20.5,200,expression(mu[0]))
text(xbar[1]+0.5,xbar[2],expression(bar(x)))

#Construct SImultan
T.ci <- function(mu, Sigma, n, avec=rep(1,length(mu)),level=0.95){
  p <- length(mu)
  if(nrow(Sigma)!=p) stop("Need length(mu) == nrow(Sigma).")
  if(ncol(Sigma)!=p) stop("Need length(mu) == ncol(Sigma).")
  if(length(avec)!=p) stop("Need length(mu) ==length(avec).")
  if(level <=0 | level >= 1) stop("Need 0 < level < 1.")
  cval <- qf(level, p, n-p) * p * (n-1) / (n-p)
  zhat <- crossprod(avec, mu)
  zvar <- crossprod(avec, Sigma %*% avec) / n
  const <- sqrt(cval * zvar)
  c(lower = zhat - const, upper = zhat + const)
} 

#Dataset Identity
S<-cov(X)
n <- nrow(X)
p <- ncol(X)

#Cek
T.ci(mu=xbar, Sigma=S, n=n, avec=c(1,0,0,0)) #CI untuk mu_0 mpg
##    lower    upper 
## 16.39689 23.78436
T.ci(mu=xbar, Sigma=S, n=n, avec=c(0,1,0,0)) #CI untuk mu_0 disp
##    lower    upper 
## 154.7637 306.6801

#Gabungan CI

#CI untuk sample besar
chi <- NULL
for(k in 1:4){
  chi <- c(chi,xbar[k] - sqrt(S[k,k]/n) * sqrt(qchisq(0.95,df=p)),
           xbar[k] + sqrt(S[k,k]/n) * sqrt(qchisq(0.95,df=p)))
}
TCI <- tCI <- bon <- NULL
alpha1<-1-0.05
alpha2<-1-(0.05/2)
alpha3 <- 1 - 0.05/(2*4)
for(k in 1:4){
avec <- rep(0, 4)
avec[k] <- 1
TCI <- c(TCI, T.ci(xbar, S, n, avec,level=alpha1))
tCI <- c(tCI, xbar[k] - sqrt(S[k,k]/n) * qt(alpha2, df=n-1),
         xbar[k] + sqrt(S[k,k]/n) * qt(alpha2, df=n-1))
bon <- c(bon,xbar[k] - sqrt(S[k,k]/n) * qt(alpha3, df=n-1),
         xbar[k] + sqrt(S[k,k]/n) * qt(alpha3, df=n-1))
chi <- c(chi,xbar[k] - sqrt(S[k,k]/n) * sqrt(qchisq(0.95,df=p)),
         xbar[k] + sqrt(S[k,k]/n) * sqrt(qchisq(0.95,df=p)))
}
rtab <- rbind(TCI, tCI, bon,chi)
round(rtab, 2)
##       mpg   mpg   disp   disp     hp     hp   wt   wt   mpg   mpg   disp   disp
## TCI 16.40 23.78 154.76 306.68 104.67 188.71 2.62 3.82 16.40 23.78 154.76 306.68
## tCI 17.92 22.26 186.04 275.41 121.97 171.41 2.86 3.57 17.92 22.26 186.04 275.41
## bon 17.27 22.92 172.62 288.82 114.55 178.83 2.76 3.68 17.27 22.92 172.62 288.82
## chi 16.81 23.37 163.24 298.21 109.35 184.02 2.68 3.75 16.81 23.37 163.24 298.21
##         hp     hp   wt   wt
## TCI 104.67 188.71 2.62 3.82
## tCI 121.97 171.41 2.86 3.57
## bon 114.55 178.83 2.76 3.68
## chi 109.35 184.02 2.68 3.75