ZADANIE 5
Kristína Mečiarová
V tejto úlohe sa pozrieme na hladiny a sily testov zo záveru 5.
cvičenia. Zvážme náhodný výber
\(\mathbf{X}_1,\dots,\mathbf{X}_n \sim
\textit{N}_2 \bigg( \begin{pmatrix} 1
\\ 2 \end{pmatrix}, \begin{pmatrix} 1 & 0.5 \\ 0.5 &
2 \end{pmatrix} \bigg)\). V našom záujme bude testovať hypotézu
\(H_0:2\mu_1−\mu_2=0\) proti
alternatíve \(H_1:2\mu_1−\mu_2 \neq0\).
Ekvivalentne môžeme použíť zápis \(H_0:\mathbb{A}\mathbf{\mu}=a\), kde \(a=0,\; \mathbb{A}= (2,-1)\) a \(\mathbf{\mu}=(1,2)^T.\) Za nulovej hypotézy
platí \[T_n=n(\mathbb{A}\overline{X}_n−a)^⊤(\mathbb{A}\Sigma\mathbb{A}^⊤)^{−1}(\mathbb{A}\overline{X}_n−a)\sim
\chi^2_1.\]
Budeme uvažovať aj prípad, že kovariančná matica \(\Sigma\) nie je známa a budeme pracovať s jej odhadom \(\mathcal{S}\). V takom prípade za platnosti nulovej hypotézy platí \[\tilde{T}_n=(n−1)(\mathbb{A}\overline{X}_n−a)^⊤(\mathbb{A}\mathcal{S}\mathbb{A}^⊤)^{−1}(\mathbb{A}\overline{X}_n−a)\sim T^2_{1,n−1}.\]
Pre overenie, že tieto testy skutočne dodržiavajú predpísanú hladinu (volíme hladinu \(\alpha=0.05\)) a sú konzistentné voči alternatíve \(H_1\), uskutočníme malú simulačnú štúdiu. Budeme voliť postupne rozsahy výberov \(n=50,100,200,500,1000\), v alternatíve budeme voliť \(a=-0,5;-0,01; 0,03;1;5\) a simulácie budeme opakovať 1000-krát.
Na prvom výstupe máme hladiny testov so známou kovariančnou maticou pre vyššie zmienené rozsahy náhodných výberov. Ako vidíme, test relatívne dobre dodržiava predpísanú hladinu aj pre nižší rozsah výberu \((n=50)\).
alpha<-0.05
library("mvtnorm")
library(greekLetters)
rozsah<-c(50,100,200,500,1000) #rozsah vyberu
opak<-1000 #pocet nezavislych opakovani
mu <- c(1, 2)
Sigma <- matrix(c(1, 0.5, 0.5, 2),2,2)
A <- c(2, -1)
###############################simulacie za H0#####################
a <- 0
#Sigma znama
vysledky.znamaH0 <- vector(mode = "list", length = 1)
set.seed(0610)
tabulka<-NULL
for (l in 1:length(rozsah)){ #rozne rozsahy
n<-rozsah[l]
hladina<-NA
zam<-numeric(opak)
for (j in 1:opak){ #opakovania testov
sample<- rmvnorm(n, mu, Sigma)
Tvalue <- n * t(A %*% apply(sample, 2, mean) - a) %*% solve(A %*% Sigma %*% A) %*% (A %*% apply(sample, 2, mean) - a)
zam[j]<-as.numeric(Tvalue>=qchisq(1-alpha,1))
}
hladina<-mean(zam)
tabulka<- cbind(tabulka,hladina)}
colnames(tabulka)<-rozsah
(vysledky.znamaH0<-tabulka)## 50 100 200 500 1000
## [1,] 0.043 0.053 0.046 0.06 0.055Na druhom výstupe máme výsledky simulácií pre testy s odhadovanou kovariančnou maticou. Hladina je opäť dobre dodržaná a výsledky sa výrazne nelíšia od výsledkov so známou kovariančnou maticou.
#Sigma odhadnuta
vysledky.neznamaH0 <- vector(mode = "list", length = 1)
set.seed(0610)
tabulka<-NULL
for (l in 1:length(rozsah)){ #rozne rozsahy
n<-rozsah[l]
hladina<-NA
zam<-numeric(opak)
for (j in 1:opak){ #opakovania testov
sample<- rmvnorm(n, mu, Sigma)
TTvalue <- (n - 1) * t(A %*% apply(sample, 2, mean) - a) %*% solve(A %*% cov(sample) %*% A) %*% (A %*% apply(sample, 2, mean) - a)
zam[j]<-as.numeric(TTvalue>=qf(1-alpha,1,n-1))
}
hladina<-mean(zam)
tabulka<- cbind(tabulka,hladina)}
colnames(tabulka)<-rozsah
(vysledky.neznamaH0<-tabulka)## 50 100 200 500 1000
## [1,] 0.042 0.053 0.052 0.06 0.055#obrazok za H0
popis<-paste("Hladiny testov")
plot(rozsah, c(vysledky.znamaH0),type="n",xlab="Rozsah n",ylab="Hladina testu",cex.lab=1.5,cex.axis=1.5,cex.main=1.5, main=popis,ylim=c(0,0.08))
axis(1,at=c(100,300,500,700,900),labels=c("100","300","500","700","900"),cex.axis=1.5)
abline(h=alpha,lty=5, col=1)
#Sigma znama
lines(rozsah,vysledky.znamaH0,lwd=2,col=3)
#Sigma odhad
lines(rozsah,vysledky.neznamaH0,lwd=2,col=4)
legend("bottomright",c(expression(paste(Sigma," známa",sep="")),expression(paste(Sigma," odhadnutá",sep=""))),bty='n',col=c(3,4),lty=c(1,1),lwd=c(2,2),cex=1.5,bg="white")
Približné dodržanie hladiny a veľmi podobné výsledky simulácií testov so známou a odhadovanou kovariančnou maticou vidíme aj na obrázku vyššie.
Ďalej sa pozrieme na silu testov pre rôzne alternatívy. Prvý výstup opäť dostávame pre simulácie so známou kovariančnou maticou a druhý pre simulácie s odhadovanou kovariančnou maticou.
###############################simulacie za H1#####################
a <- c(-0.5,-0.001,0.03,1,5)
#Sigma znama
vysledky.znamaH1 <- vector(mode = "list", length = length(a))
set.seed(0610)
tabulka<-cbind(a)
for (l in 1:length(rozsah)){ #rozne rozsahy
n<-rozsah[l]
hladina<-rep(NA,length(a))
for (i in 1:length(a)) { #rozne a
zam<-numeric(opak)
for (j in 1:opak){ #opakovania testov
sample<- rmvnorm(n, mu, Sigma)
Tvalue <- n * t(A %*% apply(sample, 2, mean) - a[i]) %*% solve(A %*% Sigma %*% A) %*% (A %*% apply(sample, 2, mean) - a[i])
zam[j]<-as.numeric(Tvalue>=qchisq(1-alpha,1))
}
hladina[i]<-mean(zam)
}
tabulka<- cbind(tabulka,hladina)
}
colnames(tabulka)<-c("a",rozsah)
(vysledky.znamaH1<-tabulka)## a 50 100 200 500 1000
## [1,] -0.500 0.458 0.706 0.937 1.000 1.000
## [2,] -0.001 0.064 0.054 0.058 0.050 0.061
## [3,] 0.030 0.055 0.062 0.049 0.069 0.077
## [4,] 1.000 0.943 0.997 1.000 1.000 1.000
## [5,] 5.000 1.000 1.000 1.000 1.000 1.000Vidíme, že so zvyšujúcim sa rozsahom výberu majú oba testy väčšiu silu, pričom so zväčšujúcim sa (v absolútnej hodnote) \(a\) naberajú oba testy tiež na sile. Opäť v prípade odhadovanej kovariančnej matice pozorujeme podobné výsleky ako v prípade známej kovariančnej matice.
#Sigma odhadnuta
vysledky.neznamaH1 <- vector(mode = "list", length = length(a))
set.seed(0610)
tabulka<-cbind(a)
for (l in 1:length(rozsah)){ #rozne rozsahy
n<-rozsah[l]
hladina<-rep(NA,length(a))
for (i in 1:length(a)) { #rozne a
zam<-numeric(opak)
for (j in 1:opak){ #opakovania testov
sample<- rmvnorm(n, mu, Sigma)
TTvalue <- (n - 1) * t(A %*% apply(sample, 2, mean) - a[i]) %*% solve(A %*% cov(sample) %*% A) %*% (A %*% apply(sample, 2, mean) - a[i])
zam[j]<-as.numeric(TTvalue>=qf(1-alpha,1,n-1))
}
hladina[i]<-mean(zam)
}
tabulka<- cbind(tabulka,hladina)
}
colnames(tabulka)<-c("a",rozsah)
(vysledky.neznamaH1<-tabulka)## a 50 100 200 500 1000
## [1,] -0.500 0.429 0.684 0.936 1.000 1.00
## [2,] -0.001 0.057 0.053 0.058 0.052 0.06
## [3,] 0.030 0.051 0.062 0.050 0.066 0.08
## [4,] 1.000 0.932 0.997 1.000 1.000 1.00
## [5,] 5.000 1.000 1.000 1.000 1.000 1.00#obrazok za H1 Sigma znama
sily<-vysledky.znamaH1[,-1]
popis<-expression(paste("Sily testov pre ",Sigma, " známu a odhadnutú",sep=""))
plot(rep(rozsah,length(a)), c(sily),type="n",xlab="Rozsah n",ylab="Sila testu",cex.lab=1.5,cex.axis=1.5,cex.main=1.5, main=popis,ylim=c(0,1.1))
axis(1,at=c(100,300,500,700,900),labels=c("100","300","500","700","900"),cex.axis=1.5)
abline(h=alpha,lty=5, col=1)
for (i in 1:length(a)){
lines(rozsah,sily[i,],lwd=2,col=i)}
sily<-vysledky.neznamaH1[,-1]
for (i in 1:length(a)){
lines(rozsah,sily[i,],lwd=2,lty=2,col=i)}
L=legend("right", legend = rep(NA,10), col=1:5, lty=c(1,1,1,1,1,2,2,2,2,2),bty='n', ncol=2, x.intersp =3.5, inset=0.02)
legend(L$rect$left-90,L$rect$top+0.08, legend = c(expression(Sigma),rep(NA,4), expression(hat(Sigma)[n]),rep(NA,4)), ncol=2, inset=0.02,bty='n',x.intersp =4)
legend(x = L$rect$left, y = L$rect$top, legend = as.character(a), col=rep(NA,5), lty=c(1,2), ncol=1, bty='n',bg = NA)
Na obrázku vyššie sú zobrazené sily testov pre rôzne uvažované \(a\), pre známu kovariančnú maticu plnou
čiarou a pre odhadovanú prerušovanou čiarou. Na obrázku pekne vidíme, že
v prípade hodnôt \(a\) blízko nule je
sila rovná hladine testu, kdežto proti alternatívam so zvyšujúcimi sa (v
absolútnej hodnote) hodnotami \(a\) sú
testy pomerne silné.
Na základe tejto malej simulačnej štúdie môžeme dospieť k záveru, že oba zmienené testy sú konzistentné proti alternatíve \(H_1:2\mu_1−\mu_2 \neq a, \; a \in \mathbb{R}\).