NMST539 | Cvičenie 5 Metóda maximálnej vierohodnosti pre mnohorozmerné normálne rozdelenie LS 2020/2021 | 30/03/21 | (online výuka) zdrojový Rmd súbor (kódovanie UTF8) Outline piateho cvičenia: teoretické opakovanie metódy maximálnej vierohodnosti; aplikácia metódy pre mnohorozmerné normálne rozdelenie; maticová algebra pre infinitesimálny počet (derivácie matíc a determinantov); teoretické príklady štatistické testy pomerom vierohodnosti; Odborná literatúra: Hardle, W. and Simar.L.: Applied Multivariate Statistical Analysis. Springer, 2015 Mardia, K., Kent, J., and Bibby, J.:Multivariate Analysis, Academic Press, 1979. Metóda maximálnej vierohodnosti je jednym zo základných nástrojov určených pre odhadovanie neznámeho parametru \(\theta \in \Omega\) (prípadne mnohorozmerného – t.j. vektorového parametru \(\boldsymbol{\theta} \in \Omega\)) v určitom predpokladanom rozdelení na základe náhodného výberu. Postup je obecne založený na tzv. vierohodnostnej funkcii: pre náhodný výber \(\boldsymbol{X}_1, \dots, \boldsymbol{X}_{n}\) z \(p \in \mathbb{N}\) rozmerného rozdelenia s hustotou \(f(\boldsymbol{x}, \boldsymbol{\theta})\) je vierohodnostná funkcia definovaná predpisom \[ L(\boldsymbol{\theta}, \mathcal{X}) = \prod_{i = 1}^n f(\boldsymbol{X}_i, \boldsymbol{\theta}), \] čo je funkcia vzľadom k argumentu \(\boldsymbol{\theta} \in \Omega\) skonštruovaná na základe náhodného výberu \(\mathcal{X} = (\boldsymbol{X}_{1}, \dots, \boldsymbol{X}_{n})^\top\). Príslušný odhad \(\widehat{\boldsymbol{\theta}}\) neznámeho vektorového parametru \(\boldsymbol{\theta}\) je definovaný ako argument, ktorý maximalizuje vierohodnostnú funkciu, t.j. hodnota \(L(\boldsymbol{\theta}, \mathcal{X})\) je maximálna možná, resp. \[ \widehat{\boldsymbol{\theta}} = \mathop{Argmax}_{\boldsymbol{\theta} \in \Omega} L(\boldsymbol{\theta}, \mathcal{X}). \] Maximálne vierohodný odhad \(\widehat{\boldsymbol{\theta}}\) neznámeho vektorového parametru \(\boldsymbol{\theta} \in \Omega\) je za určitých predpokladov regularity asymptoticky normálny a platí, že \[ \sqrt{n} (\widehat{\boldsymbol{\theta}} - \boldsymbol{\theta}) \sim N_{p}(\boldsymbol{0}, \mathbb{I}_{p}^{-1}(\boldsymbol{\theta})), \quad \textrm{pre} \quad n \to \infty \] kde \(\mathbb{I}_{p}(\boldsymbol{\theta})\) je tzv. Fisherova informačná matica o neznámom vektorovom parametri \(\boldsymbol{\theta} \in \Omega\). Maximálne vierohodný odhad \(\widehat{\boldsymbol{\theta}}\) je teda asymptotický nestranný a eficientný (v tom zmysle, že asymptotická variančná-kovariančná matica \(\mathbb{I}_{p}^{-1}(\boldsymbol{\theta})\) dosahuje Rao-Cramérovu úroveň). 1. Maticová algebra infinitesimálneho počtu V následujúcej časti stručne pripomenieme niektoré užitočné pravidla pre maticové počítane v súvislosti s deriváciou matíc a determinantov. V prípade mnohorozmerného parametru \(\boldsymbol{\theta} = (\theta_1, \dots, \theta_p)^\top \in \Omega \subseteq \mathbb{R}^p\) je výhodné pracovať s pomocou tohto maticového kalkulu. Pomocou matíc menších rozmerov (napr. \(2 \times 2\), alebo \(3 \times 3\)) samostatne overte platnosť nasledujúcich pravidiel: Príklad 1: Pre obecnú štvorcovú maticu \(\mathcal{X}\) typu \(p \times p\) platí: \[\frac{\partial |\mathcal{X}|}{\partial \mathcal{X}} = \big(Adj(\mathcal{X})\big)^\top;\] Príklad 2: Pre symetrickú štvorcovú maticu \(\mathcal{X}\) typu \(p \times p\) platí: \[\frac{\partial |\mathcal{X}|}{\partial \mathcal{X}} = 2 (Adj(\mathcal{X})\big) - Diag(Adj(\mathcal{X}));\] Príklad 3: Pre obecné štvorcové matice \(\mathcal{X}\) a \(\mathbb{A}\) typu \(p \times p\) platí: \[\frac{\partial tr(\mathcal{X} \mathbb{A})}{\partial \mathcal{X}} = \frac{\partial tr(\mathbb{A} \mathcal{X})}{\partial \mathcal{X}} = \mathbb{A}^\top;\] Príklad 4: Pre štvorcové matice \(\mathcal{X}\) a \(\mathbb{A}\) typu \(p \times p\), kde \(\mathbb{X}\) je navyše symetrická, platí \[\frac{\partial tr(\mathcal{X} \mathbb{A})}{\partial \mathcal{X}} = \frac{\partial tr(\mathbb{A} \mathcal{X})}{\partial \mathcal{X}} = 2 \mathbb{A} - Diag(\mathbb{A});\] Príklad 5: Pre obecné štvorcové matice \(\mathcal{X}\) a \(\mathbb{A}\) typu \(p \times p\), platí: \[\frac{\partial tr(\mathcal{X}^{-1} \mathbb{A})}{\partial \mathcal{X}} = \frac{\partial tr(\mathbb{A} \mathcal{X}^{-1})}{\partial \mathcal{X}} = - \mathcal{X}^{-1}\mathbb{A}^\top \mathcal{X}^{-1};\] Príklad 6: Pre matice \(\mathbb{A}\) a \(\mathbb{B}\) typu \(q \times p\) a \(p \times q\) a štvorcovú maticu \(\mathcal{X}\) typu \(p \times p\), platí: \[\frac{\partial tr(\mathbb{A}\mathcal{X} \mathbb{B})}{\partial \mathcal{X}} = \frac{\partial tr(\mathbb{B} \mathbb{A} \mathcal{X})}{\partial \mathcal{X}} = \mathbb{B}\mathbb{A};\] Príklad 7: Pre obecnú štvorcovú maticu \(\mathcal{X}\) typu \(p \times p\) platí: \[\frac{\partial \ln |\mathcal{X}|}{\partial \mathcal{X}} = \frac{1}{det(\mathcal{X})} Adj(\mathcal{X}) = (\mathcal{X}^{-1})^\top \equiv \mathcal{X}^{{-\top}};\] Príklad 8: Pre symetrickú štvorcovú maticu \(\mathcal{X}\) typu \(p \times p\) platí: \[\frac{\partial \ln |\mathcal{X}|}{\partial \mathcal{X}} = 2 \mathcal{X}^{-1} - Diag(\mathcal{X}^{-1});\] Uvedené pravidlá sú užitočné pri derivácii vierohodnostnej funkcie a získanie príslušného maximálne vierohodného odhadu v rôznych mnohorozmerných rozdeleniach. V následujúcej časti niektoré z týchto pravidiel využijeme pri konštrukcii maximálne vierohodného odhadu neznámeho vektorového parametru pre rôzne rozdelenia. 2. Metóda maximálnej vierohodnosti pre mnohorozmerný parameter Pre prvé dva príklady je pre získanie maximálne vierohodného odhadu postačujúce využiť štandardný postup na základe parciálných derivácii. Pre posledné dva príklady je ale omnoho výhodnejšie využiť vhodné vzťahy uvedené vyššie. Príklad 9 Uvažujme dvojrozmerný náhodný výber \(\boldsymbol{X}_1, \dots, \boldsymbol{X}_{n}\), kde \(\boldsymbol{X}_{i} = (X_{i 1}, X_{i, 2})^\top\), z rozdelenia s hustotou \[ f_{\boldsymbol{\theta}}(x_{1}, x_{2}) = \frac{1}{\theta_1\theta_2} \cdot exp\Big\{ -\big( \frac{x_1}{\theta_1} + \frac{x_2}{\theta_2} \big) \Big\} \cdot \mathbb{I}_{\{x_1 > 0\}} \mathbb{I}_{\{x_2 > 0\}}, \] kde \(\boldsymbol{\theta} = (\theta_1, \theta_2)^\top \in (0, \infty) \times (0, \infty)\). Nájdite odhad neznámeho vektorového parametru pomocou metódy maximálnej vierohodnosti a spočítajte Rao-Cramérovu mez. Príklad 10 Uvažujte náhodný výber z rozdelenia, ktoré je definované dvojromernou hustotou \[ f_{\boldsymbol{\theta}}(x_1, x_2) = \frac{1}{\theta_1^2 \theta_2} \frac{1}{x_2} exp\Big\{ - \big( \frac{x_1}{\theta_1 x_2} + \frac{x_2}{\theta_1 \theta_2} \big) \Big\}, \] pre \(x_1, x_2 > 0\) a hustota je nulová inak. Nájdite maximálne vierohodný odhad neznámeho vektorového parametru \(\boldsymbol{\theta} = (\theta_1, \theta_2)^\top \in (0, \infty) \times (0, \infty)\) a spočítajte Rao-Cramérovu mez pre \(\widehat{\boldsymbol{\theta}}\). Príklad 11 Uvažujme náhodný vektor \(\boldsymbol{X}\) s \(p\) rozmerným normálnym rozdelením \(N_{p}(\boldsymbol{\mu}, \Sigma)\), kde \(\boldsymbol{\mu} \in \mathbb{R}^p\) a \(\Sigma = Diag\{\sigma_{11}, \dots, \sigma_{pp}\}\), pre \(\sigma_{jj} > 0\) a všetky \(j = 1, \dots, p\). Nech \(\boldsymbol{X}_{1}, \dots, \boldsymbol{X}_{n}\) predstavuje náhodný výber z daného rozdelenia. Pomocou náhodného výberu a metódy maximálnej vierohodnosti nájdite odhad pre neznámy vektor stredných hodnôt \(\boldsymbol{\mu}\) a pre variančnú-kovariančnú maticu \(\Sigma\). Odvoďte Rao-Cramérovu mez pre neznámy vektorový parameter \(\boldsymbol{\theta} = (\mu_1, \dots, \mu_p, \sigma_{11}, \dots, \sigma_{pp})^\top\). Príklad 12 Uvažujte mnohorozmerný náhodný výber \(\boldsymbol{X}_1, \dots, \boldsymbol{X}_n\) z mnohorozmerného normálneho rozdelenia \(N_{p}(\boldsymbol{\mu}, \Sigma)\) s neznámym vektorom stredných hodnôt \(\mu \in \boldsymbol{\mu} \in \mathbb{R}^p\) a neznámou symetrickou a pozitívne definitnou variančnou-kovariačnou maticou \(\Sigma \in \mathbb{R}^{p \times p}\). Pomocou metódy maximálnej vierohodnosti nájdite odhad pre neznámy vektorový parameter \(\boldsymbol{\mu}\) a variančnu-kovariančnu maticu \(\Sigma\). 3. Štatistické testy pomerom vierohodnosti Example 13 Uvažujte náhodný výber \(\boldsymbol{X}_1, \dots, \boldsymbol{X}_n\) z mnohorozmerného normálneho rozdelenia \(N_{p}(\boldsymbol{\mu}, \Sigma)\), kde \(\boldsymbol{\mu} \in \mathbb{R}^p\) je neznámy vektor stredných hodnôt \(\Sigma\) je neznáma variančná-kovariančná matica. Pre konkrétnu voľbu \(\boldsymbol{\mu}_0 \in \mathbb{R}^p\) nás zaujíma test nulovej hypotézy \[ H_{0}: \boldsymbol{\mu} = \boldsymbol{\mu}_{0}, \] oproti obecnej (obojstrannej) alternatíve \[ H_{1}: \boldsymbol{\mu} \neq \boldsymbol{\mu}_{0}. \] Ukážte, že za platnosti nulovej hypotézy je maximum vierohodnostnej funkcie rovné \(\ell_{0} = \ell(\mathcal{X}, \mu_0, \mathcal{S} + (\overline{\boldsymbol{X}}_n - \boldsymbol{\mu}_0)(\overline{\boldsymbol{X}}_{n} - \boldsymbol{\mu}_0)^\top))\). Ukážte, že testová štatistika založená na pomere vierohodnosti je \[ - 2 \log \lambda = 2 (\ell_1 - \ell_0) = n \log\Big(1 + (\overline{\boldsymbol{X}}_n - \boldsymbol{\mu}_{0}) ^\top \mathcal{S}^{-1} (\overline{\boldsymbol{X}}_n - \boldsymbol{\mu}_{0})\Big). \] Testová štatistika závisí na \((\overline{\boldsymbol{X}}_n - \boldsymbol{\mu}_{0})^\top \mathcal{S}^{-1} (\overline{\boldsymbol{X}}_n - \boldsymbol{\mu}_{0})\) a navyše, za platnosti nulovej hypotézy \(H_0\) platí, že \[ (n - 1)(\overline{\boldsymbol{X}}_n - \boldsymbol{\mu}_{0}) ^\top \mathcal{S}^{-1} (\overline{\boldsymbol{X}}_n - \boldsymbol{\mu}_{0}) \sim T_{p, n - 1}^2. \] S využitím znalosti o Hotellingovom \(T^2\) rozdelení a jeho súvislosti s Fisherovým F rozdelením (viď minulé čvičenie) je vpodstate už triviálne daný test v programe R naprogramovať a aplikovať na konkrétne data. V nasledujúcom príklade sa zameriame na obecnejšiu variantu testu. Príklad 14 V prvom kroku budeme simulovať náhodný výber z konkrétneho dvojrozmerného normálneho rozdelenia \(N_{2}\Big(\Big(\begin{array}{c}1\\2\end{array}\Big), \Big( \begin{array}{cc}1 & 0.5\\0.5 & 2 \end{array} \Big)\Big)\). Následne, pomocou vygenerovaného náhodného výberu, chceme testovať nulovú hypotézu (oproti obecnej alternatíve) \[ H_{0}: 2\mu_{1} - \mu_2 = 0.\] Zaujíma nás teda konkrétna lineárna kombinácia jednotlivých zložiek neznámeho vektoru stredných hodnôt. Rozlíšime dva prípady: variančná-kovariančná matica \(\Sigma\) je známa, prípadne \(\Sigma\) je neznáma matica. library("mvtnorm") n <- 100 mu <- c(1, 2) Sigma <- matrix(c(1, 0.5, 0.5, 2),2,2) set.seed(1234) ### simulacia nahodneho vyberu sample <- rmvnorm(n, mu, Sigma) V obecnom tvare je možné nulovú hypotézu vyjadriť v tvare \[ H_0: \mathbb{A}\boldsymbol{\mu} = a, \] kde \(\boldsymbol{\mu} = (\mu_1, \mu_2)^\top = (1, 2)^\top\), \(\mathbb{A} = (2, -1)\), a \(a = 0\). Obecná alternatíva má tvar \[ H_1: \mathbb{A}\boldsymbol{\mu} \neq a, \] V prípade, že variančná-kovariančná matica \(\Sigma\) je známa, ľahko sa overí, že za platnosti nulovej hypotézy dostaneme \[ T_{n} = n(\mathbb{A}\overline{\boldsymbol{X}}_{n} - a)^\top(\mathbb{A}\Sigma \mathbb{A}^\top)^{-1} (\mathbb{A} \overline{\boldsymbol{X}}_{n} - a) \sim \chi_{1}^{2}. \] a nulová hypotéza je zamietaná pre veľké hodnoty testovej štatistiky \(T_{n}\) (t.j., pre hodnoty \(T_{n}\), ktoré sú väčšie, než príslušný \((1 - \alpha)\) kvantil \(\chi_2\) rozdelenia – v našom konkrétnom prípade s jedným stupňom voľnosti, kde \(\alpha \in (0, 1)\)). Pomocou nasledujúcich príkazov spočítame hodnotu testovej štatistiky pre nasimulované data: a <- 0 A <- c(2, -1) (Tvalue <- n * t(A %*% apply(sample, 2, mean) - a) %*% solve(A %*% Sigma %*% A) %*% (A %*% apply(sample, 2, mean) - a)) ## [,1] ## [1,] 0.3746285 a túto hodnotu je potrebné porovnať s príslušným kvantilom \(\chi_1^2(0.95) = 3.8415\). Nulovú hypotézu teda na hladine \(\alpha = 0.05\) nezamietame (keďže hodnota testovej štatistiky \(T_n\) je približne \(0.37\), čo je menej, ako daná kritická hodnota). Ak budeme ale predpokladať, že variančná-kovariačná matica \(\Sigma\) je tiež neznáma, musíme test založiť na inej testovej štatistike, ktorú jednoducho získame tak, že neznámu maticu nahradíme príslušným (empirickým) odhadom. Samotné rozdelenie testovej štatistiky za platnosti nulovej hypotézy sa ale zmení, keďže tentokrat je nutné počítať aj s určítou mierou variability, ktorá je prítomná v odhade neznámej variančnej-kovariančnej matice. Opäť priamočiaro dostaneme, že platí \[ \tilde{T}_{n} = (n - 1)(\mathbb{A}\overline{\boldsymbol{X}}_{n} - a)^\top(\mathbb{A}\mathcal{S} \mathbb{A}^\top)^{-1} (\mathbb{A} \overline{\boldsymbol{X}}_{n} - a) \sim T_{1, n - 1}^2, \] pričom konkrétnu hodnotu spočítame pomocou programu R (TTvalue <- (n - 1) * t(A %*% apply(sample, 2, mean) - a) %*% solve(A %*% cov(sample) %*% A) %*% (A %*% apply(sample, 2, mean) - a)) ## [,1] ## [1,] 0.3735325 a túto hodnotu tentokrát porovnáme s príslušným kvantilom Fisherovho rozdelenia \(F_{1, 99}(0.95) = 3.9371\). Nulovú hypotézu teda na rovnakej hladine \(\alpha = 0.05\) zamietame. Oba testy vychádzaju podobne – a to hlavne v prípadoch, keď je rozsah náhodného výberu dostatočne veľký a odhad variančnej-kovariančnej matice je konzistentný. Domáca (samostatná) úloha (Deadline: Šieste cvičenie / 06.04.2021) Rozmyslite si malú simulačnú štúdiu, v ktorej sa zameriate na empirické vlastnosti štatistického testu pomerom vierohodnosti pri známej, resp. neznámej variančnej-kovariačnej matici v mnohorozmernom normálnom rozdeleni (alternatívne môžete uvažovať aj test pomerom vierohodnosti v nejakom inom rozdeleni). Data simulujte za platnosti nulovej hypotézy a pre oba prípady (t.j., známa vs. neznáma variančná-kovariančná matica) sledujte pravdepodobnosť chyby prvého druhu (t.j. počet prípadov, kedy bola nulová hypotéza zamietnutá voči celkovému počtu uskutočnených testov). Generujte data za platnosti vhodnej (t.j. konkrétnej) alternatívnej hypotézy a sledujte sílu testu. Využijte rôzne nastavenia pre veľkosť náhodného výberu, špecifikáciu alternatívnej hypotézy, prípadne aj variančnú-kovariančnú maticu. Riešenie umiestnite na svoju webovú stránku, najneskôr v utorok, 06.04.2021, do 14:00.