Testy o střední hodnotě při známém rozptylu

Mgr. Zdeňka Geršlová

Teorie

Prezentace Statistical inference I and II, kapitola 7.

Příklady

Příklad 1

Vzájemné porovnání testovacích statistik \(U_W, U_S, U_{LR}\)

Vzájemně porovnejte tvary křivek \(y = \frac{x}{1 + x}\), \(y = \log(1 + x)\) a \(y = x\) a stanovte, která křivka dosahuje na intervalu \(\left(0, 5 \right)\) (a) nejvyšších hodnot (b) nejnižších hodnot. Výsledek porovnání lze využít při porovnávání testovacích statistik.

Postup: Vytvoříme sekvenci bodů na zadaném intervalu, definujeme jednotlivé křivky a vykreslíme je do grafu pomocí funkce lines.

Výsledek

Figure 1: Porovnání křivek (globální a lokální pohled)

Příklad 2

Aktuální vs. nominální hladina významnosti \(\alpha\), konzervativní vs. liberální test

Nechť

  1. \(X \sim N(20, 100)\)
  2. \(X \sim p N(20, 100) + (1 - p) N(28, 100)\), kde \(p = 0.9\), tedy jde o směs dvou normálních rozdělení \(X \sim N(20, 100)\) a \(X \sim N(28, 100)\) v poměru \(9:1\).

Pro obě části (a) i (b) vygenerujte \(M = 500\) náhodných výběrů s rozsahem \(n = 5\), resp. \(n = 100\) a vypočítejte hodnotu testovací statistiky \(Z_W\) pro Waldův test nulové hypotézy \(H_0: \mu = \mu_0 = 20\) oproti \(H_1: \mu \neq 20\), když \(\sigma ^ 2\) známe (\(\sigma ^ 2 = 10 ^ 2\)) na hladině významnosti \(\alpha = 0.05\). Hodnoty testovacích statistik \(Z_W\) zaneste do histogramu. Vždy spočítejte, kolik testovacích statistik \(Z_W\) spadá do kritického oboru \(W\). Podíl tohoto čísla a \(M\) představuje aktuální hladinu významnosti \(\widehat{\alpha}\). Porovnejte ji s nominální hladinou významnosti \(\alpha\) a v každé ze čtyř situací určete, zda je test konzervativní či liberální.

Pozn.: Test samozřejmě nemusí být ani konzervativní ani liberální, což je ten nejlepší případ. Konzervativnost nebo liberálnost testu chápeme jako (negativní) vlastnost testu, kterou, je-li přítomná, musíme mít na zřeteli.

Postup

Pro univerzální použití vytvoříme funkci HistWald, která bude připravena pro libovolnou směs normálních rozdělení s volitelnými parametry jednotlivých rozdělení, počtu a rozsahu náh. výběrů.

HistWald <- function(mu0, n, M = 500,
                   mu1 = 20, mu2 = mu1, 
                   sigma1 = 10, sigma2 = sigma1, 
                   p = 0.9, alpha = 0.05, main) {
  X <- matrix(NA, M, n)
for (i in 1:M) {
  bin <- rbinom(n, 1, p)
  X[i, ][bin == 1] <- rnorm(sum(bin), mu1, sigma1)
  X[i, ][bin == 0] <- rnorm(n - sum(bin), mu2, sigma2)
} # generovani smesi rozdeleni
  m <- apply(X, 1, mean)
  zW <- ... # vzorec pro Zw statistiku
  d <- hist(zW, plot = F)$dens
  xfit <- seq(min(zW) - 10, max(zW) + 10, length = 512) # sekvence pro vykresleni Zw 
  yfit <- ... # hustota N(0,1)
  alpha_aktual <- sum(abs(zW) > qnorm(1 - alpha / 2)) / M # aktualni hl. vyzn.
  d <- hist(zW, plot = F)$dens # vyska sloupcu histogramu
hist(..., ylim = c(0, max(yfit, d)), ...)
mtext(expression(z[W]), side = 1, line = 2.2)
mtext(bquote(paste(alpha == .(alpha), "; ", widehat(alpha) == .(alpha_aktual))), side = 1,
      line = 3.3)
mtext(main, side = 1, line = 4.5)
lines(...) # krivka hustoty 
}

Výsledek pro normální rozdělení

source("M8986-source.R")
set.seed(10)
par(mar = c(6, 4, 1, 1), mfrow = c(1,2))
HistWald(mu0 = 20, n = 5,  main = expression(paste('X ~ N(20, 100); n = 5'))) 
HistWald(mu0 = 20, n = 100, main = expression(paste('X ~ N(20, 100); n = 100')))

Rozdělení Waldovy testovací statistiky pro test o střední hodnotě při známém rozptylu

Směs normálních rozdělení

par(mar = c(6, 4, 1, 1), mfrow = c(1,2))
HistWald(mu0 = 20, mu2 = 28, n = 5,
         main = expression(paste('X ~ 0.9 N(20, 100) + 0.1 N(28, 100); n = 5')))

HistWald(mu0 = 20, n = 100, mu2 = 28, 
       main = expression(paste('X ~ 0.9 N(20, 100) + 0.1 N(28, 100); n = 100')))

Rozdělení Waldovy testovací statistiky pro test o střední hodnotě při známém rozptylu (směs rozdělení)

Závěr

Při opakovaném spuštění MC studie vidíme, že v případě (a) aktuální hladina významnosti \(\widehat{\alpha}\) nabývá rovnoměrně častokrát vyšší i nižší hodnoty než je nominální hladina významnosti \(\alpha\) (navíc rozdíly nejsou příliš velké). DIS pro \(\mu\), když \(\sigma ^ 2\) známe, není v tomto případě konzervativní ani liberální, stejně jako test založený na testovací statistice \(Z_W\).

Oproti tomu v případě (b) je aktuální hladina významnosti \(\widehat{\alpha}\) vyšší než nominální hladina významnosti \(\alpha\), tj. aktuální pravděpodobnost pokrytí je nižší než nominální pravděpodobnost pokrytí a DIS pro \(\mu\), když \(\sigma ^ 2\) známe, je v tomto případě liberální, stejně jako test založený na testovací statistice \(Z_W\). Efekt se více projeví pro vyšší rozsah náhodných výběrů.

Příklad 3

Rozdělení testovací statistiky pro test o střední hodnotě se známým rozptylem

Nechť náhodný výběr \(X\) pochází z normálního rozdělení, t.j. \(X \sim N(\mu, \sigma ^ 2)\). Pomocí simulační studie porovnejte rozdělení testovací statistiky \(Z_W\) pro test nulové hypotézy \(H_0\): \(\mu = 150\) (alternativní hypotéza \(H_1\): \(\mu \neq 150)\), když rozptyl \(\sigma ^ 2\) známe, s rozdělením testovací statistiky stanovené na základě náhodného výběru se střední hodnotou \(\mu\). Parametry zvolte
(a) \(\mu = 146\), \(\sigma ^ 2 = 10 ^ 2\), \(n = 50\);
(b) \(\mu = 155\), \(\sigma ^ 2 = 10 ^ 2\), \(n = 50\).

Nechť dále \(X\) pochází ze směsi dvou normálních rozdělení, t.j. \(X \sim [p N(\mu, 10 ^ 2) + (1 - p) N(\mu, 30 ^ 2)]\), kde \(p = 0.9\) a
(c) \(\mu = 146\);
(d) \(\mu = 155\).
Proveďte simulační studii popsanou výše také pro tento náhodný výběr.

Postup

  1. Nasimulujte M pseudonáhodných výběrů, \(M = 1, \dots, 2\,000\) a pro každý vypočítejte realizaci testovací statistiky \(z_{W, \lambda} ^ {(m)}=\frac{\bar{x}_m - \mu_0}{\sigma} \sqrt{n}\) pro nulovou hypotézu
    \(H_0\): \(\mu = 150\) oproti \(H_1\): \(\mu \neq 150\).
  2. Vykreslete histogram testovacích statistik \(Z_W\) a superponujte jej jednak křivkou hustoty normálního rozdělení \(N(\lambda, 1)\) s parametrem necentrality \(\lambda\) \((\lambda = \frac{\mu - \mu_0}{\sigma} \sqrt{n}\), kde \(\mu\) je skutečná střední hodnota (relevantní za platnosti \(H_1\))) a jednak křivkou hustoty standardizovaného normálního rozdělení \(N(0, 1)\).

Obě křivky nakonec vzájemně porovnejte.

Pozn.: U směsi křivka hustoty necentrálního rozdělení nesuperponuje histogram statistik \(z_W\) dostatečně. Zamyslete se nad tím, proč.

Výsledné grafy

par(mar = c(6, 4, 1, 1), mfrow = c(2,2)) 
RozdeleniNecentr(mu = 146, main = 'X ~ N(146, 100)')
RozdeleniNecentr(mu = 155, main = 'X ~ N(155, 100)', pozice = 'topright') 
RozdeleniNecentr(mu = 146, sigma2 = 30, main = 'X ~ 0.9 N(146, 100) + 0.1 N(146, 900)') 
RozdeleniNecentr(mu = 155, sigma2 = 30, main = 'X ~ 0.9 N(155, 100) + 0.1 N(155, 900)',
                 pozice = 'topright')

V source vytvoříme funkci RozdeleniNecentr, která bude pro směs rozdělení vykreslovat histogramy superponované křivkami hustoty (normálního rozdělení s parametrem necentrality a standardizovaného norm. rozdělení).

RozdeleniNecentr <- function(mu, mu0 = 150, sigma = 10,
                             sigma2 = sigma, M = 2000, n = 50,
                             p = 0.9, main = "", pozice = "topleft"){
X <- ... # matice nah. vyberu - analogicky predchozimu pr.
m <- apply(X, 1, mean)
zW <- ... # test. statistika z_w
lambda <- ... # parametr necentrality

xfit <- ... # sekvence pro vykresleni hustoty
yfit <- ... # hustota norm. rozd.
zfit <- ... # hustota necentralniho norm. rozd.

hist(...)
lines(...) # cervene N(0,1), modre necentralni
legend(pozice, ...)
}

Komentář

Centrální rozdělení přísluší parametru \(\mu_0 = 150\), proto se červená křivka centrálního rozdělení realizuje okolo hodnoty \(0\). Naproti tomu rozdělení testovací statistiky \(Z_W\) přísluší hodnotě \(\mu = 146\) (resp. \(\mu = 155\)), v grafu znázorněno jako modrý histogram superponovaný modrou křivkou, proto parametr necentrality \(\lambda\) nabývá záporné (resp. kladné) hodnoty a histogram i s modrou křivkou necentrálního rozdělení se realizuje nalevo (resp. napravo) od červené křivky centrálního rozdělení. Modrá křivka superponuje histogram přesně, protože náhodný výběr pochází z předpokládaného rozdělení \(N(146, 10 ^ 2)\).

V případě směsi rozdělení je situace podobná, ovšem modrá křivka nyní nesuperponuje histogram přesně, což má důvod právě v přítomnosti “příměsi” v námi předpokládaném rozdělení.

Příklad 4

Silofunkce pro jednovýběrový Z-test o střední hodnotě

Předpokládejme, že \(X \sim N(\mu, \sigma ^ 2)\), kde \(\sigma ^ 2\) známe. Nechť \(\theta = \mu\). Na hladině významnosti \(\alpha = 0.05\) testujeme všechny tři typy hypotéz

  1. \(H_{01}:\) \(\mu = \mu_0\) oproti \(H_{11}:\) \(\mu \neq \mu_0\) (oboustranná);
  2. \(H_{02}:\) \(\mu \leq \mu_0\) oproti \(H_{12:}\) \(\mu > \mu_0\) (pravostranná);
  3. \(H_{03}:\) \(\mu \geq \mu_0\) oproti \(H_{13}:\) \(\mu < \mu_0\) (levostranná).

Odvoďte tvary silofunkcí pro všechny tři typy hypotéz (a)–(c), t.j. tvary \(\beta_{11} ^ * (\mu)\), \(\beta_{12} ^ * (\mu)\) a \(\beta_{13} ^ * (\mu)\) (viz Aplikovaná štatistická inferencia I - př. 155, str. 122 a dále)

Dále nakreslete silofunkce pro všechny tři typy hypotéz (a)–(c), kde \(\mu_0 = 150\), a \(\sigma ^ 2 = 10 ^ 2\). Do jednoho obrázku zakreslete vždy tvary silofunkcí pro \(n = 10\), \(n = 20\), \(n = 50\) a \(n = 100\). Hodnoty \(\mu\) volte rozumně, např. v intervalu \(\left( 132; 168\right)\).

Postup

Vytvoříme funkci SilaExakt pro výpočet exaktní silofunkce

SilaExakt <- function(mu0 = 0, mu, sigma = 1,
                      n, alpha = 0.05, alternative = "two.sided") {
  if (alternative == "two.sided") {
    sila <- pnorm(qnorm(alpha / 2) - (mu - mu0) / sigma * sqrt(n)) +
      pnorm(qnorm(alpha / 2) + (mu - mu0) / sigma * sqrt(n))
  }
  if (alternative == "greater") {   ... # pravostranna alternativa
    }
  if (alternative == "less") {  ... # levostranna alternativa
  }
  return(sila)
}

A dále funkci PlotSila, která bude pracovat s výstupem této funkce a vykreslí silofunkce pro zadaná n.

PlotSila <- function(mu0 = 150, mu, sigma = 10, 
                     n, alpha = 0.05, alternative = "two.sided") {
  barva <- ... # zde muzeme definovat posloupnost barev
  plot(mu, SilaExakt(...), type = "n", ylim = c(0, 1), ...)

  for (i in 1:length(n)) {
    sila <- SilaExakt(..., n = n[i],...)
    lines(mu, sila, col = barva[i])
  }
  legend(..., col = barva, legend = paste("n =", n), ... ) # legenda
  abline(...) # pridani linie v alpha
}
mu <- seq(..., length = 512) # sekvence pro mu
n <- c(...) # vektor n

Výsledek

Silofunkce testu o střední hodnotě se známým rozptylem (oboustranná, pravostranná a levostranná alternativa)

Komentář

Hladina významnosti \(\alpha\) (z definice) udává riziko chyby, že \(H_0\) nesprávně zamítáme, přestože platí. Ve všech třech případech (a), (b) i (c) se tedy pro \(\mu = \mu_0 = 150\) hodnota síly rovná přímo hodnotě hladiny významnosti \(\alpha = 0.05\).

(a) Z grafu pozorujeme, jak s rostoucí vzdáleností \(\mu\) od \(\mu_0=150\) roste pravděpodobnost, že \(H_0\) zamítáme (tj. roste síla testu). Platí tedy, že síla testu klesá s hodnotou \(\mu \rightarrow \mu_0 = 150\) z obou stran.

(b) U pravostranné alternativy síla testu klesá s \(\mu \rightarrow \mu_0 = 150\) zprava.

(c) U levostranné alternativy síla testu klesá s hodnotou \(\mu \rightarrow \mu_0 = 150\) zleva.

Příklad 5

Porovnání exaktní a aproximatické silofunkce
Uveďte tvary přesné silofunkce \(\beta_{11} ^ *\) a přibližné silofunkce \(\tilde{\beta}_{11} ^ *\) pro test \(H_{01}:\) \(\mu = \mu_0\) oproti \(H_{11}:\) \(\mu \neq \mu_0\) když \(\sigma ^ 2\) známe. Nakreslete křivky obou silofunkcí do jednoho grafu, kde na ose \(x\) budou různé hodnoty parametru \(\mu\) na ose \(y\) vynesená silofunkce, a porovnejte jejich tvary. Výsledek slovně okomentujte.

Hodnotu \(n\) zvolte 100, \(\mu_0 = 150\) a \(\sigma ^ 2 = 10 ^ 2\). Rozsah osy \(x\) volte rozumně, pro globální pohled např. \(\langle 145; 155\rangle\), pro lokální zaměření rozdílů zvolte rozsah osy \(x\) \(\langle 148; 152 \rangle\).

Pozn.: Aproximatickou sílu můžeme spočítat pouze pro oboustrannou alternativu, neboť její myšlenka spočívá ve společném vyjádření obou částí síly (dvou distribučních funkcí) prostřednictvím jedné distribuční funkce s absolutní hodnotou. U jednostranných alternativ je síla tvořena pouze jednou distribuční funkcí, proto zde myšlenka fungující u oboustranné alternativy postrádá smysl.

Vytvoříme funkci pro výpočet aproximatické síly SilaAprox a pro vykreslení použijeme plot() a lines()

SilaAprox <- function(...){
  sila <- pnorm(qnorm(alpha / 2) + abs(mu0 - mu) / sigma * sqrt(n))
}

# nastaveni pro lokalni pohled
plot(..., xlim = c(148, 152), ylim = c(-0.1, 0.4), asp = F, las = 1) 

Výsledné grafy

Exaktní vs. aproximatická silofunkce testu o střední hodnotě při známém rozptylu

Závěr

Exaktní a empirická silofunkce jsou si tvarově velmi blízké všude s výjimkou intervalu cca \((149.2 ; 150.8)\). Exaktní sílu tedy můžeme aproximovat v případě, že \(\mu\) je dostatečně vzdálená od \(\mu_0\).

Důvod: Aproximatická síla je založena na zanedbání jedné ze dvou distribučních funkcí, které tvoří exaktní sílu. V okolí \(\mu = \mu_0 = 150\) je vzdálenost \(\mu\) od \(\mu_0\) malá a do výsledné hodnoty exaktní síly přispívají velkým dílem obě distribuční funkce. Pokud tedy jednu z těchto distribučních funkcí zanedbáváme, přicházíme v aproximatické síle o její příspěvek (aproximatická síla je v okolí \(\mu = \mu_0\) výrazněji menší než exaktní síla).

Příklad 6

Minimální rozsah náhodného výběru
Předpokládejme, že \(X \sim N(\mu, \sigma ^ 2)\), kde \(\sigma ^ 2 = 10 ^ 2\). Nechť \(\theta = \mu\). Testujeme všechny tři typy hypotéz:

\(H_{01}:\) \(\mu = \mu_0\) oproti \(H_{11}:\) \(\mu \neq \mu_0\) (oboustranná),
\(H_{02}:\) \(\mu \leq \mu_0\) oproti \(H_{12:}\) \(\mu > \mu_0\) (pravostranná),
\(H_{03}:\) \(\mu \geq \mu_0\) oproti \(H_{13}:\) \(\mu < \mu_0\) (levostranná),

kde \(\mu_0 = 150\). Vypočítejte minimální rozsah náhodného výběru pro test nulové hypotézy při \(\alpha = 0.05\) a \(1 - \beta = 0.8\), pro
\(\mu \in \{145, 145.5, \dots 154.5, 155\}\) (ad (1)); \(\mu \in \{ 150.25, 150.5, 150.75, \dots 153.5, 153.75, 154 \}\) (ad (2)); \(\mu \in \{ 146, 146.25, 146.5,\) \(\dots 149.5, 149.75 \}\) (ad (3)).

  • Závislost minimálního rozsahu náhodného výběru na hodnotě \(\mu\) zakreslete do grafu pomocí bodů (na osu \(x\) vyneste parametr \(\mu\), na osu \(y\) minimální rozsah náhodného výběru).

  • Sestavte tabulku minimálních rozsahů náhodného výběru pro test nulové hypotézy \(H_0: \mu = \mu_0\), kde \(\mu_0 = 150\) oproti alternativním hypotézám \(H_{11}\), \(H_{12}\) a \(H_{13}\) při předem stanovené síle \(\beta ^ * = 0.8\) a hladině významnosti \(\alpha = 0.05\), předpokládáme-li, že výběrová střední hodnota \(\mu\) bude nabývat hodnot \(\mu \in \{145, 146, 147, 148, 149, 149.5, 150.5, 151,\$ 152, 153, 155\}\).

Postup v R

Sestavíme funkci MinZTest, která bude pro zadané parametry rozdělení a stanovenou sílu a hladinu významnosti počítat minimální rozsah pro jednotlivé alternativy (nezapomeňte výstup zaokrouhlit na nejbližší vyšší celé číslo).

MinZTest <- function(mu, mu0 = 150, sigma = 10,
                      alpha = 0.05, sila = 0.8, alternative = "two.sided") {
  if (alternative == "two.sided") {
    n <- (qnorm(sila) - qnorm(alpha / 2))^2 / abs(mu - mu0)^2 * sigma^2
  }
  ... # jednostranne alternativy
  return(ceiling(n)) # zaokrouhleni
}

Lze použít i funkci power.z.test z knihovny asbio.

library(asbio)
result <- power.z.test(sigma = 10, power = 0.8, alpha = 0.05, 
                       test = "two.tail", effect = c(0.5, 1, 2))
ceiling(result$n)
[1] 3140  785  197

Potom sestavíme sekvenci \(\mu\) dle zadání pro příslušnou alternativu a vykreslíme bodový graf, kde na ose \(x\) jsou jednotlivé hodnoty \(\mu\) a na ose \(y\) příslušný minimální rozsah.

Pozn.: K vykreslení čárkované čáry v místě, kde nemá smysl uvažovat min. rozsah (případ jednostranné alternativy), lze použít funkci segments.

segments(145, 0 , 150 , 0, lty = 2, col = 'grey40')

Grafy

Minimální rozsahy náh. výběru pro Z test o střední hodnotě při známém rozptylu

Tabulka

```{r}
#| code-line-numbers: false
#| tbl-cap: Rozsahy náh. výběru pro jednotlivé alternativy Z testu pro zvolená $\\mu$
#| options: knitr.kable.NA = '.'
library(knitr)
mu_vyber <- c(145:149, 149.5, 150.5, 151:155)
n11 <- MinZTest(mu_vyber, alternative = 'two.sided')
n12 <- MinZTest(mu_vyber, alternative = 'greater')
n13 <- MinZTest(mu_vyber, alternative = 'less')
n12 <- c(rep(NA, 6), n12[7:12])
n13 <- c(n13[1:6], rep(NA, 6))

tab <- data.frame(t(data.frame(mu = mu_vyber, n11 = n11, n12 = n12, n13 = n13)), 
                    row.names = c('$\\mu$', '$H_{11}: \\mu = \\mu_0$',
                                  '$H_{12}: \\mu > \\mu_0$', '$H_{13}: \\mu < \\mu_0$'))
kable(tab, col.names = NULL)
```
Rozsahy náh. výběru pro jednotlivé alternativy Z testu pro zvolená \(\mu\)
\(\mu\) 145 146 147 148 149 149.5 150.5 151 152 153 154 155
\(H_{11}: \mu = \mu_0\) 32 50 88 197 785 3140.0 3140.0 785 197 88 50 32
\(H_{12}: \mu > \mu_0\) NA NA NA NA NA NA 2474.0 619 155 69 39 25
\(H_{13}: \mu < \mu_0\) 25 39 69 155 619 2474.0 NA NA NA NA NA NA

Pozn.: Pro oboustrannou alternativu jsou rozsahy náhodných výběrů (pro pevně zvolené \(\mu_0\) a \(\mu\)) vyšší než pro jednostranné alternativy.

Příklad 7

Silofunkce testu o střední hodnotě \(\mu\) když \(\sigma ^ 2\) známe
Předpokládejme, že \(X \sim N(\mu, \sigma ^ 2)\), kde \(\sigma ^ 2 = 10 ^ 2\), \(n = 100\). Nechť \(\theta = \mu\). Na hladině významnosti \(\alpha = 0.05\) testujeme hypotézu \(H_{01}: \mu = \mu_0\) oproti \(H_{11}: \mu \neq \mu_0\) (oboustranná), kde \(\mu_0 = 150\).

Vytvořte animaci zobrazující
(a) změnu polohy necentrálního rozdělení vzhledem k hodnotě centrálního rozdělení testovací statistiky testu o \(\mu\) když \(\sigma ^ 2\) známe, spolu s barevně odlišenou oblastí kritického oboru a exaktní síly;
(b) změnu hodnoty exaktní silofunkce; při měnící se střední hodnotě náhodného výběru \(\mu = 140, 141,\dots, 146, 146.5, \dots,\) \(153.5, 154, 155, \dots 160\).

Pozn: V prvním grafu budeme mít červeně křivku hustoty centrálního rozdělení s vyznačeným kritickým oborem, modře křivku hustoty necentrálního rozdělení s parametrem necentrality lambda a vyznačenou exaktní silou testu (síla - pravděpodobnost - obsah pod křivkou vymezený kvantilem).

Animace

Průběh síly testu pro střední hodnotu při známém rozptylu, obooustranná alternativa

Animace - kód

```{r}
#| fig-show: animate
#| animation-hook: gifski
#| fig-cap: Průběh síly testu pro střední hodnotu při známém rozptylu,
#|  obooustranná alternativa
#| code-line-numbers: false
mu <- c(140:145, seq(146, 153.5, by = 0.5), 154:160) 
for (i in 1:length(mu)) {
  SilaAnimace(150, mu = mu[i], sigma = 10, n = 100,
              alternative = 'two.sided')
}
```

Průběh síly testu pro střední hodnotu při známém rozptylu, obooustranná alternativa

Funkce SilaAnimace

Vytvoříme funkci SilaAnimace, která bude zobrazovat oba grafy.
Nejprve vypočítáme parametr necentrality lambda, určíme dostatečně hustou posloupnost x, ve které vykreslujeme hustotu, a vypočítáme hustotu centrálního a necentrálního rozdělení.

SilaAnimace <- function(mu0, mu, sigma, n,
                        alpha = 0.05,
                        alternative = 'two.sided') {
  lambda <- (mu - mu0) / sigma * sqrt(n)
  x <- ... # posloupnost na intervalu (-10,10)
  y <- ... # hustota centralniho rozdeleni
  l <- ...(..., mean = lambda) # hustota necentralniho rozd.
  
  ## kod pro vykresleni grafu
}

Vykreslení hustoty

Pro vykreslení hustoty použijeme nejprve plot s argumentem type = 'n' pro nastavení grafického okna a potom postupně vykreslíme křivky hustot pomocí lines a oblast kritického oboru pomocí funkce polygon (na vstupu je nejprve vektor x-ových souřadnic a potom y-ových souřadnic, níže je postup pro křivku hustoty centrálního rozdělení, necentrální zobrazíme analogicky, jen použijeme hustotu necentrálního rozdělení a zobrazíme modrou barvou).

par(mfrow = c(1, 2), mar = c(5, 4, 1, 1))
plot(..., type = 'n', las = 1)

q1 <- qnorm(alpha / 2)
q2 <- qnorm(1 - alpha / 2)
  
polygon(x = c(min(x), max(x[x <= q1]),
          max(x[x <= q1]), rev(x[x <= q1])), 
          y = c(0, 0, y[x == max(x[x <= q1])],
          rev(y[x %in% x[x <= q1]])), 
          col = 'red', density = 40)
polygon(c(min(x[x >= q2]), max(x),
          rev(x[x >= q2]), min(x[x >= q2])),
          c(0, 0, rev(y[x %in% x[x >= q2]]),
          y[x == min(x[x >= q2])]),
          col = 'red', density = 40)
lines(x, y, col = 'red')

Polygon

Vykreslení síly

Ve druhém grafu vykreslíme křivku exaktní síly (pomocí funkce SilaExakt, kterou jsme vytvořili v předchozím příkladu) pro rozumnou sekvenci středních hodnot a přidáme vždy výrazný bod v místě aktuální hodnoty \(\mu\). Vodorovnou čarou zvýrazníme hodnotu hladiny významnosti \(\alpha\).

mu1 <- seq(mu0 - 20, mu0 + 20, length = 512)

sila_ex <- SilaExakt(mu0 = mu0, mu = mu1, sigma = sigma,
                       n = n, alternative = 'two.sided')
plot(mu1, sila_ex, ... )

sila_akt <- SilaExakt(mu0 = mu0, mu = mu, sigma = sigma,
                      n = n, alternative = 'two.sided')
points(...)

Závěr

V hodnotě \(\mu = \mu_0 = 150\) necentrální rozdělení splývá s centrálním a silofunkce dosahuje svého minima, které je rovno hladině významnosti \(\alpha = 0.05\).

Při vzdalování \(\mu\) od \(\mu_0\) se necentrální rozdělení vzdaluje od centrálního a hodnota silofunkce roste.