DÚ 5

V tejto úlohe budeme generovať náhodný výber z rozdelenia \(\mathbf{X}_1,\dots,\mathbf{X}_n \sim \mathrm{N}_2 \bigg( \begin{pmatrix} 2 \\ 2 \end{pmatrix}, \begin{pmatrix} 1 & 0.5 \\ 0.5 & 2 \end{pmatrix} \bigg)\).

library("mvtnorm")

n <- 100
mu <-  c(2, 2)
Sigma <- matrix(c(1, 0.5, 0.5, 2),2,2)

set.seed(1111)
sample <- rmvnorm(n, mu, Sigma)

Budeme testovať hypotézu \(H_0: \mu_1 - \mu_2 = 0\) voči alternatíve \(H_1: \mu_1 - \mu_2 \neq 0.\) Tento zápis vieme ekvivalentne napísať ako \(H_0:\mathbb{A}\boldsymbol{\mu}=a\), kde \(a=0,\; \mathbb{A}= (1,-1)\) a \(\boldsymbol{\mu}=(2,2)^\top\). Pre testovanie rozlíšime dva prípady. Najprv budeme predpokladať, že matica \(\Sigma\) je známa a potom budeme pracovať s jej odhadom \(\mathcal{S}\).

Za platnosti nulovej hypotézy a známej kovariančnej matice dostávame testovú štatistiku

\(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.\)

Nulovú hypotézu zamietame ak platí \(T_{n} \geq \chi_1^2(0.95) = 3.842\).

A <- c(1, -1) 
(Tvalue <- n * t(A %*% apply(sample, 2, mean)) %*% solve(A %*% Sigma %*% A) %*% (A %*% apply(sample, 2, mean)))
##           [,1]
## [1,] 0.2117112

Hodnota testovej štatistiky je 0.212 a teda nezamietame nulovú hypotézu v prospech alternatívy.

Ďalej predpokladajme, že maticu \(\Sigma\) nepoznáme a budeme pracovať s jej odhadom \(\mathcal{S}\). Testová štatistika má za platnosti nulovej hypotézy tvar

\(\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}\).

Tentokrát porovnávame testové štatistiku s príslušným kvantilom Fisherovho rozdelenia t. j. zamietame nulovú hypotézu v prospech alternatívy, ak \(\tilde{T}_{n} \geq F_{1, 99}(0.95) = 3.937\).

(Tvalue2 <- (n - 1) * t(A %*% apply(sample, 2, mean)) %*% solve(A %*% cov(sample) %*% A) %*% (A %*% apply(sample, 2, mean)))
##           [,1]
## [1,] 0.1633993

Hodnota testovej štatistiky je 0.163 a opäť nezamietame nulovú hypotézu.

Ďalej spočítame pravdepodobnosť chyby prvého druhu (pravdepodobnosť zamietnutia platnej hypotézy) pre rôzne veľkosti náhodného výberu. Budeme pracovať s výbermi veľkosti \(n=50, 100, 500, 1000, 10\,000\) a simulácie budeme opakovať 1000 krát.

##          n   pst
## [1,]    50 0.046
## [2,]   100 0.056
## [3,]   500 0.041
## [4,]  1000 0.053
## [5,] 10000 0.046

Z vyššie uvedenej tabuľky vidíme, že testy dobre dodržiavajú predpísanú hladinu. Pozrime sa, či dostaneme takisto dobré výsledky aj za predpokladu, že maticu \(\Sigma\) nepoznáme.

##          n  pst2
## [1,]    50 0.046
## [2,]   100 0.046
## [3,]   500 0.049
## [4,]  1000 0.050
## [5,] 10000 0.041
x=c(1:5)
plot(x,pst,xaxt="n",xlab="n=50,100,500,1 000,10 000",col="red",pch=16,ylim=c(0,0.1),ylab="Hladina testu",main="Chyba prvého druhu")
lines(x,pst,col="red",type = "b")
lines(x, pst2, pch = 16, col = "blue",type = "b")
abline(h=0.05,col="black",lty=2)

legend("topleft", legend=c("známa kovariančná matica", "odhadnutá kovariančná matica"),
       col=c("red", "blue"), lty = 1:1, cex=0.8)

Vidíme, že výsledky sú v poriadku aj v prípade, keď musíme počitať s odhadom \(\mathcal{S}\) matice \(\Sigma\).

Ostáva spočítať pravdepodobnosť chyby druhého druhu (pravdepodobnosť nezamietnutia neplatnej hypotézy) pre rôzne alternatívy aj rozsahy výberov. Voľme preto \(a = -0.01, 0.5, 1\) a rovnako ako doteraz, \(n=50, 100, 500, 1000, 10\,000\). Opäť spočítame tieto pravdepodobnosti pre známu kovariančnú maticu ako aj pre jej odhad.

##          n a=-0.01 a=0.5   a=1
## [1,]    50   0.043 0.718 0.999
## [2,]   100   0.058 0.933 1.000
## [3,]   500   0.049 1.000 1.000
## [4,]  1000   0.055 1.000 1.000
## [5,] 10000   0.134 1.000 1.000

Podobne pre odhadnutú kovariančnú maticu \(\mathcal{S}\).

##          n a=-0.01 a=0.5   a=1
## [1,]    50   0.044 0.704 0.999
## [2,]   100   0.038 0.949 1.000
## [3,]   500   0.049 1.000 1.000
## [4,]  1000   0.060 1.000 1.000
## [5,] 10000   0.100 1.000 1.000
x=c(1:5)
plot(x,pst[,1],xaxt="n",xlab="n=50,100,500,1 000,10 000",col="red",pch=16,ylim=c(0,1.1),ylab="Sila testu",main="Chyba druhého druhu")
lines(x,pst[,1],col="red",type = "b")
lines(x,pst[,2],col="red",type = "b", lty = 4)
lines(x,pst[,3],col="red",type = "b", lty = 5)

lines(x, pst2[,1], pch = 16, col = "blue",type = "b")
lines(x, pst2[,2], pch = 16, col = "blue",type = "b", lty = 4)
lines(x, pst2[,3], pch = 16, col = "blue",type = "b", lty = 5)

abline(h=1,col="black",lty=2)

legend("right", legend=c("známa kovariančná matica, a=-0.01","známa kovariančná matica, a=0.5","známa kovariančná matica, a=1", "odhadnutá kovariančná matica, a=-0.01","odhadnutá kovariančná matica, a=0.5","odhadnutá kovariančná matica, a=1"),
       col=c("red","red","red","blue","blue","blue"), cex=0.8, lty = c(1,4,5,1,4,5))

Vidíme, že pre známu aj odhadnutú kovariančnú maticu je sila testu pre \(a = 0.5\) aj \(a = 1\) veľmi dobrá.