Ústav matematiky a statistiky Přírodovědecká fakulta Masarykova univerzita Štatistická inferencia I Zadania príkladov a domácich úloh a niektoré riešenia Stanislav Katina katina@math.muni.cz 14. decembra 2014 Obsah 1 Model rozdelenia pravdepodobnosti a štatistický model 1 2 Charakteristiky polohy a variability a štatistická grafika 19 3 Testovanie hypotéz 20 Katina, S., 2014: Štatistická inferencia I 1 1 Model rozdelenia pravdepodobnosti a štatistický model Príklad 1 (porovnanie dvoch typov modelov) Model rozdelenia pravdepodobností je modelom náhodnej premennej X, napr. model rozdelenia pravdepodobnosti náhodnej premennej X šírka dolnej čeľuste alebo (2) model rozdelenia pravdepodobnosti náhodnej premennej X hrúbka kožných rias u dospelých zdravých žien. Statistický model je modelom náhodnej premennej Y\X (Y kauzálne závisí na X), napr. (1) model závislosti náhodnej premennej Y šírka dolnej čeľuste na premennej X pohlavie alebo (2) model závislosti náhodnej premennej Y hrúbka kožných rias u dospelých zdravých žien na premennej X BMI. Všimnime si, že náhodné premenné označujeme X alebo Y podľa toho, aký model ich charakterizuje. pred Príklad 2 (jednoduchý náhodný výber) V jednoduchom náhodnom výbere s rozsahom n z populácie s konečným rozsahom N má každý prvok rovnakú pravdepodobnosť vybratia. Ak vyberáme bez vrátenia, hovoríme o jednoduchom náhodnom výbere bez vrátenia1 (Dalgaard 2008). Ak vyberáme s vrátením, hovoríme o jednoduchom náhodnom výbere s vrátením2. Majme množinu M. s N = 10 prvkami a chceme z nej vybrať n = 3 prvkov (a) bez vrátenia a (b) s vrátením. Koľko máme možností? Ako vyzerá jedna takáto možnosť, ak ide o množinu Ai = {1, 2,..., 10}. Zopakujte to isté pre N = 100, n = 30 a množinu .M = {1, 2,..., 100}. cvič Riešenie aj v Q§ (a) Spolu máme (^) možných náhodných výberov. Ak N = 10 a n = 3, potom kombinačné číslo (n) = (NNn)M = (3°) = 120 možností. Ak N = 100 a n = 30, potom (^) = f™) = 2.937234 x 1025 možností. choose(10,3) # počet všetkých možných výberov bez vrátenia choose(100,30) library(utils) combn(10,3) # počet všetkých možných výberov bez vrátenia combn(100,30) sample(x=l:10,size=3,replace = FALSE) # jednoduchý náhodný vyber bez vrátenia sample(x=l:100,size=30,replace = FALSE) (b) Spolu máme ^N+n~ľ^ možných náhodných výberov. Ak N = 10 a n = 3, potom ^N+n~ľ^ = ^0 = (10+3-1) = 220 možností. Ak N = 100 a n = 30, potom f^1) = (100+f-1) = 2.009491 x 1029 možností. choose(10+3-1,3) # počet všetkých možných výberov s vratenim choose(100+30-1,30) library(utils) combn(10 + 3-1 ,3) # počet všetkých možných výberov s vratenim combn(100+30-1,30) sample(x=l:10,size=3,replace = TRUE) # jednoduchý náhodný vyber s vratenim sample(x=l:100,size=30,replace = TRUE) Príklad 3 (jednoduchý náhodný výber) Nech je skupina ľudí označená identifikačnými číslami (ID) od 1 do 30. Vyberte (a) náhodne 5 ľudí z 30 bez návratu, (b) náhodne 5 ľudí z 30 s návratom a nakoniec (c) náhodne 5 ľudí z 30 bez návratu, kde ľudia s ID od 28 do 30 majú pravdepodobnosť vybratia 4x väčšiu ako ľudia s ID od 1 do 27. cvič. 1 Kombinácie bez opakovania n-tej triedy z N prvkov množiny A4. 2Kombinácie s opakovaním n-tej triedy z N prvkov množiny A4. (14. decembra 2014) Katina, S., 2014: Štatistická inferencia I 2 Riešenie v <® 1 sample(x=l:30,size=5,replace = FALSE) 2 sample(x=l:30,size=5,replace = TRUE) 3 sample(x=l:30,size=5, prob=c(rep(1/39,27),rep(4/39,3)), replace = FALSE) Príklad 4 (normálne rozdelenie) Majme náhodnú premennú X (môže to byť napr. výška postavy 10-ročných dievčat) a predpokladáme, že má normálne rozdelenie s parametrami fi (stredná hodnota) a a2 (rozptyl), čo zapisujeme ako X ~ N(fi,a2), fi = 140.83, a2 = 33.79. Normálne rozdelenie predstavuje model rozdelenia pravdepodobnosti pre túto náhodnú premennú. Vypočítajte pravdepodobnosť Pr (a < X < b) = Pr (X < b)-Pr(X < a) = Fx (b)-Fx (a), kde a = fi-ka, b = /j + ka, k = 1,2,3.3 Nakreslite hustotu rozdelenia pravdepodobnosti, vyfarbite oblasť medzi bodmi a a b a popíšte osi x a y tak, ako je uvedené na obrázku 1. cvič Riešenie (aj v <3i); (pozri obrázok 1) a = fi-a = 135.0171, b = fi + a = 146.6429, Pr (\X - fi\ > a) = 0.3173, Pr (\X - fi\ < a) = 1 - 0.3173 = 0.6827, a = /i-2a= 129.2042, b = /i + 2a = 152.4558, Pr (\X - fi\ > 2a) = 0.0455, Pr (\X - fj,\ < 2a) = 1 - 0.0455 = 0.9545, a = p - 3a = 123.3913, b = + 3a = 158.2687, Pr (\X - /x| > 3a) = 0.0027, Pr (\X - /x| < 3a) = 1 - 0.0027 = 0.9973. Alternatívny výpočet cez štandardizované normálne rozdelenie (syn. normálne normované rozdelenie) je nasledovný: mu <- 0 sig <- 1 bin <- seq(mu-3 *sig,mu + 3*sig,by = sig) pnorm(bin[7]) - pnorm(bin [1]) # 0.9973002 pnorm(bin[6]) - pnorm(bin [2]) # 0.9544997 pnorm(bin[5]) - pnorm(bin [3]) # 0.6826895 Dostaneme pravidlo 68.27 — 95.45 — 99.73 (tzv. „miery normálneho rozdelenia"). Príklad 5 (normálne rozdelenie) Majme X ~ N(fi,a2), kde fi = 150, a2 = 6.25. Vypočítajte a = fi — Xi_aa a b = fi + X\_aa tak, aby Pr (a < X < b) = 1 — a, bola rovná 0.90, 0.95 a 0.99. Číslo xi-a Je kvantil normálneho normovaného rozdelenia, t.j. Pľ(Z = < Xi_a) = 1 — a, Z ~ N(0,1). Nakreslite hustotu rozdelenia pravdepodobnosti, vyfarbite oblasť medzi bodmi a a b a popíšte osi x a y tak, ako je uvedené na obrázku 2. cvič Dostaneme pravidlo 90 — 95 — 99 (tzv. „upravené miery normálneho rozdelenia"). Použili sme nerovnosť Pr (ua/2 < Z < u1_a/2) = $ {u\-a/2) — $ (ua/2) = 1 — a, kde $ je distribučná funkcia normálneho normovaného rozdelenia a všeobecne a G (0,1/2); v príklade a = 0.1, 0.05 a 0.01. Pozri obrázok 2. Príklad 6 (normálne rozdelenie) Predpokladajme model normálneho rozdelenia N (132,132) pre systolický krvný tlak. Aká časť populácie (v %) bude mať hodnoty väčšie ako 160 mm Hg? cvic 3Pravdepodobnosť Pr (a < X < b) — Pr (a < X < b), pretože pravdepobnosť v bode (tu a a b) je rovná nule pre spojité premenné, t.j. Pr(a) — Pr(6) — 0. Pre diskrétne premenné to neplatí. (14. decembra 2014) Katina, S., 2014: Štatistická inferencia I 3 normovaná výska normovaná výska normovaná výska Pr(-1.6449 16) = l-£ťMi<15Pr (X = Xi) = l-£ťMi<15 ©^'(l-p)""*' = l"IU 16), ale my potrebujeme Pr(X > 16). Preto 16) = 1 - E^x ^ ^- pred Príklad 13 (štandardizované normálne rozdelenie) Model pre náhodný výber X1} X2,..., Xn je ÍV(0,1) a hovoríme, že X±,X2,... ,Xn pochádza zo štandardizovaného normálneho rozdelenia, t.j. X ~ N(pi,a2), kde pi = 0 a a2 = 1. Parameter modelu N(pi,a2) je vektor 0 = (0,1). Hustota tohto rozdelenia má tvar 0 (x) = f (x) = -^ž^e 2,i:Gl. pred Príklad 14 (dvojrozmerné normálne rozdelenie) Náhodný vektor (X, Y)T má dvojrozmerné normálne rozdelenie N2 (fi, S) , kde [i = (pLľ,pL2)T a S = 4Aproximácia znamená „približné vyjadrenie", t.j. buď nejaké rozdelenie aproximujeme iným (majúcim isté výhody oproti tomu, ktoré aproximujeme), alebo aproximujeme dáta nejakým rozdelením (ktoré popisuje dáta pomocou l'ahko interpretovatelhých parametrov). / a\ paľa2\ \paia2 o\ J (14. decembra 2014) Katina, S., 2014: Štatistická inferencia I 5 i-1-1-1-1-r i-1-1-1-1-r Obr. 3: Aproximácia binomického rozdelenia normálnym pre p = 0.515 a N = 5,10 a 50; spojnicový graf superponovaný hustotou (prvý riadok) a distribučnou funkciou (druhý riadok) s hustotou f (x, y) exp 2(1-P2) „2 2p (71(72 „2 kde (x,y) G K2, Pj G K1, o"2 > 0, j = 1,2, p G (—1,1) sú parametre, potom 6 = (//i, pi2, p)-Výraz v exponente môžeme písať ako x - jjli \ í a{ p&i&2 \ I x - p,i y - P2 J °l J \y-p2 marginálne rozdelenia5 sú X ~ N (p±,af) aY ~ N (p2, 0"|)> P Je koeficient korelácie® (pozri obrázok 5). cvič. Príklad 15 (dvojrozmerné normálne rozdelenie) (1) Nakreslite hustotu dvojrozmerného normálneho rozdelenia N2 (^t, S) pomocou funkcie image() a superponujte ho s konturovým grafom hustoty toho istého rozdelenia pomocou funkcie contourO. (2) Nakreslite hustotu dvojrozmerného normálneho rozdelenia N2 (n, S) pomocou funkcie persp(). Hustotu rozsekajte na 12 intervalov, kde hodnoty v týchto intervaloch budú zodpovedať farbám terrain. colors (12). Použite nasledovné parametre (a) fi1 = 0,p2 = 0, o x = l,a2 = 1, p = 0; (b) px = 0,p2 = 0, (Ji = l,o-2 = 1, P = 0.5; (c) pi = 0,p2 = 0, (7i = 1, a2 = 1.2, p = 0.5. Vzorové riešenie pozri na obrázku 5. cvič. 5Margináne rozdelenie je rozdelenie marginálnej náhodnej premennej, tu X nezávisle na y a naopak Y nezávisle na X. 6Z tohto príkladu je zrejmé, že na dostatočný popis dvojrozmerného normálneho rozdelenia potrebujeme päť parametrov, t.j. strednú hodnotu a rozptyl pre marginálne rozdelenie náhodných premenných X s. Y s. korelačný koeficient p — p(X,Y) popisujúci silu lineárneho vzťahu X a Y. (14. decembra 2014) Katina, S., 2014: Štatistická inferencia I 6 Bin(5,0.1] Bin(10,0.i; 1-1-1 0 5 10 Bin(50,0.1] Bin(5,0.i; 5 Bin(10,0.i; —I— ro Bin(50,0.1] Obr. 4: Aproximácia binomického rozdelenia normálnym pre p = 0.1 a N = 5,10 a 50; spojnicový graf superponovaný hustotou (prvý riadok) a distribučnou funkciou (druhý riadok) Príklad 16 (dvojrozmerné normálne rozdelenie) Nech náhodnou premennou X je najväčšia výška mozgovne u mužov (skull.pH; v mm) a náhodnou premennou Y je morfologická výška tváre u mužov (face.H; v mm); dáta: one-sample-correlation-skull. txt. Nech E [X] = p\ je stredná hodnota najväčšej výšky mozgovne a Var[X] = o\ je rozptyl najväčšej výšky mozgovne, E[Y] = p2 je stredná hodnota morfologickej výšky tváre a Var[Y] = a2 je rozptyl morfologickej výšky tváre. Predpokladajme, že najväčšia výška mozgovne X má normálne rozdelenie N(pi,o-\) a morfologická výška tváre Y má normálne rozdelenie N(p2,o-^). Potom (X, Y)T má dvojrozmerné normálne rozdelenie N2 (^t, S) s parametrami [i = (pi, p2)T, čo je vektor stredných hodnôt a o\, o\ a p, čo sú parametre kovariančnej matice S, kde sila lineárneho vzťahu týchto dvoch premenných je daná veľkosťou a znamienkom p. Potom 0 = (pi,p2,o~f,p)T■ (1) Nakreslite hustotu dvojrozmerného normálneho rozdelenia N2 (n, S) pomocou funkcie image() a superponujte ho s konturovým grafom hustoty toho istého rozdelenia pomocou funkcie contourO. (2) Nakreslite dvojrozmerný jadrový odhad hustoty pomocou funkcií kde2d() a image() a superponujte ho s konturovým grafom hustoty dvojrozmerného normálneho rozdelenia N2 (fi, S) pomocou funkcie contourO. Hustotu rozsekajte na 12 intervalov, kde hodnoty v týchto intervaloch budú zodpovedať farbám terrain. colors (12). Namiesto 0 použite vektor 0 odhadnutý z dát. Riešenie pozri na obrázku 6. cvič. Príklad 17 (štandardizované dvojrozmerné normálne rozdelenie) Náhodný vektor {X, Y)T má dvojrozmerné normálne rozdelenie N2 (0,S), kdeO = (0,0)T a S 1P s hustotou 1 í x2 - 2pxy + y2 0 ("'V) = f ^ V) = ^T=T2 \ " 2(1-p>) (14. decembra 2014) Katina, S., 2014: Štatistická inferencia I 7 |ii = 0, |i2 = O, CT-i = 1, a2 = 1, p = O |i! = O, |i2 = O, CT-i = 1, a2 = 1, p = 0.5 |ii = O, |i2 = O, ctí ctí o 'en o .p CO _ i i i i r 122 132 142 najväčšia výska mozgovne (mm) ctí > ctí ctí o 'en o .p í I I I I I 122 132 142 najväčšia výska mozgovne (mm) Obr. 6: Hustota dvojrozmerného normálneho rozdelenia s parametrom 0, ktorý je odhadnutý z dát (vľavo) a superimpozícia kontúr hustoty dvojrozmerného normálneho rozdelenia s parametrom 0. ktorý je odhadnutý z dát a dvojrozmerného jadrového odhadu hustoty (vpravo) 2) knižnice library(mvtnorm) a funkcie rmvnormO; 3) funkcie rnorm() a nasledovného algoritmu - nech X\ ~ N(0,1) a X2 ~ N(0,1); potom (Y1,Y2) ~ N2 (n, S), kde [i = (//i, p2)T, čo je vektor stredných hodnôt a o\, a2 a p, čo sú parametre kovariančnej matice S, kde sila lineárneho vzťahu Y\ a Y2 je daná veíkosíou a znamienkom p; Y\ = o\X\ + p\ a Y2 = a2(pXi + a/1 — p2X2) + p2. Nasimulujte pseudonáhodné čísla Y\ a Y2 z N2 (fi, S). Vypočítajte dvojrozmerný jadrový odhad hustoty (Y"1; Y"2)T pomocou funkcie kde2d(). Nakreslite ho pomocou funkcie imageO a superponujte ho s konturovým grafom hustoty dvojrozmerného normálneho rozdelenia N2 (fi, S) pomocou funkcie contourO. Pri simulácii použite nasledovné parametre (a) p1 = 0,p2 = 0, (Ji = 1,(72 = 1, p = 0; (1) n = 50 a (2) n = 1000; (b) px = 0,p2 = 0, (7i = 1,(72 = 1, p = 0.5; (1) n = 50 a (2) n = 1000; (c) pľ = 0,p2 = 0, (7i = 1,(72 = 1-2, p = 0.5; (1) n = 50 a (2) n = 1000. Vzorové riešenie pozri na obrázku 6. Príklad 20 (dvojrozmerné normálne rozdelenie) Simuláciu pseudonáhodných čísel z N2 (fi, S) môžeme v Oř urobiť použitím nasledovných alternatívnych funkcií: 1) knižnice library(mass) a funkcie mvrnormO; 2) knižnice library(mvtnorm) a funkcie rmvnormO; 3) funkcie rnorm() a nasledovného algoritmu - nech X\ ~ N(0,1) a X2 ~ ÍV(0,1); potom (Yi,Y2) ~ N2 (fi, S); kde [i = (pi, p2)T, čo je vektor stredných hodnôt a o\, a2 a p, čo sú parametre kovariančnej matice S; kde sila lineárneho vzťahu Y\ a Y2 je daná veľkosťou a znamienkom p; Y\ = o\X\ + p\ a Y2 = a2(pXi + a/1 — p2X2) + p2. Nasimulujte pseudonáhodné čísla Y\ a Y2 z N2 (fi, S). Vypočítajte dvojrozmerný jadrový odhad hustoty {Yi,Y2)T pomocou funkcie kde2d(). Nakreslite ho pomocou funkcie imageO a superponujte ho s konturovým grafom hustoty dvojrozmerného normálneho rozdelenia N2 (fi, S) pomocou funkcie contourO. Hustotu rozsekajte na 12 intervalov, kde hodnoty v týchto intervaloch budú zodpovedať farbám terrain. colors (12). Pri simulácii použite nasledovné parametre (a) p1 = 0,p2 = 0, (7i = 1,(72 = 1, p = 0; (1) n = 50 a (2) n = 1000; (b) p1 = 0,p2 = 0, (7i = 1,(72 = 1, p = 0.5; (1) n = 50 a (2) n = 1000; (c) pľ = 0,p2 = 0, (7i = 1,(72 = 1-2, p = 0.5; (1) n = 50 a (2) n = 1000. Vzorové riešenie pozri na obrázku 7. (14. decembra 2014) Katina, S., 2014: Štatistická inferencia I 9 Obr. 7: Hustoty dvojrozmerného normálneho rozdelenia (prvý riadok n = 50; druhý riadok n = 1000) Príklad 21 (zmes dvoch dvojrozmerných normálnych rozdelení) Simuláciu p s eudonáhodný ch čísel zo zmesi dvoch normálnych rozdelení pN2 (^1;Si) + (1 — p)N2 (fi2,'E2) môžeme v ® urobit použitím jedného z alternatínych postupov z príkladu 20. Nasimulujte pseudonáhodné čísla X a Y (í) zo zmesi pN2 + (1 -p)N2 {il2,H2), kde 0 = (pn, /i12, o\x, af2, pu /i21, p22, o~L P2Y « (2) z dvojrozmerného rozdelenia N2 (^t, S), kde parametre predstavujú spoločný vektor stredných hodnôt a spoločnú kovariančnú maticu, t.j. 0 = (pi, P2, o"2, o"!; P)T■ ^re W vypočítajte dvojrozmerný jadrový odhad hustoty (X, Y)T pomocou funkcie kde2d(). (a) Nakreslite teoretickú hustotu (2) pomocou funkcie image() a superponujte ju s konturovým grafom teoretickej hustoty (2) pomocou funkcie contourO. (b) Nakreslite teoretickú hustotu (1) pomocou funkcie image() a superponujte ju s konturovým grafom teoretickej hustoty (1) pomocou funkcie contourO■ (c) Nakreslite dvojrozmerný jadrový odhad hustoty realizací (1) pomocou funkcie image() a superponujte ju s konturovým grafom teoretickej hustoty (1) pomocou funkcie contourO■ Hustotu rozsekajte na 12 intervalov, kde hodnoty v týchto intervaloch budú zodpovedať farbám terra-in. colors (12). Pri simulácii použite 0 = (—1.2, —1.2,1,1, 0,1,1,1,1, 0)T, (1) 0 = (/in, /Í12, & 11, & 12, Pi, P21, P22, ofi) of2, P2)T; ni = n>2 = 50 a p = 0.5 (odhady pochádzajú z nasimulovaných dát). (2) 0 = (pi, P2, cr2, a2, p)T a n\ = n2 = 50 (odhady pochádzajú zo spoločného výberu nasimulovaných dát). Vzorové riešenie pozri na obrázku 8. DU Príklad 22 (zmes dvoch dvojrozmerných normálnych rozdelení) Nech (Xi,Yí)T pochádza z rozdlenia N2 (fi1, Si), kde X\ je priemerná dĺžka dolnej končatiny lowex.L v milimetroch a Y\ dĺžka trupu tru.L v milimetroch (u mužov). Nech (X2,l2)T pochádza z rozdlenia N2 (fi2,'Ľ2), kde X2 je (14. decembra 2014) Katina, S., 2014: Štatistická inferencia I 10 --1-1-1-1-1-1-1-1 f H-1-1-1-1-1-1-1-1 f H-1-1-1-1-1-1-1-1 -4-3-2-1012 3 4 -4 -3 -2 -1 0 1 2 3 4 -4 -3 -2 -1 0 1 2 3 4 Obr. 8: Spoločná hustota dvojrozmerného normálneho rozdelenia (vľavo), hustota zmesi dvoch dvojrozmerných normálnych rozdelení (uprostred) a dvojrozmerný jadrový odhad superponovaný hustotou zmesi dvoch dvojrozmerných normálnych rozdelení (vpravo) priemerná dĺžka dolnej končatiny lowex.Ĺ v milimetroch a Y2 dĺžka trupu tru.Ĺ v milimetroch (u žien). Predpokladajme, že X je priemerná dĺžka dolnej končatiny a Y dĺžka trupu pochádzajú (1) zo zmesi pN2(fi1,'S1) + (1 - p)N2 (fi2, £2), kde 0 = (fiu, fi12, a^, af2, pľ, fi21, fi22, a%2, p2)T a (2) z dvojrozmerného rozdelenia N2 (fi, S), kde parametre predstavujú spoločnú vektor stredných hodnôt a spoločnú kovariančnú maticu, t.j. 0 = /x2, o"2, a2, p)T. Pre (1) vypočítajte dvojrozmerný jadrový odhad hustoty (X, Y)T pomocou funkcie kde2d(). (a) Nakreslite teoretickú hustotu (2) pomocou funkcie image() a superponujte ju s konturovým grafom teoretickej hustoty (2) pomocou funkcie contourO. (b) Nakreslite teoretickú hustotu (1) pomocou funkcie image() a superponujte ju s konturovým grafom teoretickej hustoty (1) pomocou funkcie contourO■ (c) Nakreslite dvojrozmerný jadrový odhad hustoty realizací (1) pomocou funkcie image() a superponujte ju s konturovým grafom teoretickej hustoty (1) pomocou funkcie contourO■ Hustotu rozsekajte na 12 intervalov, kde hodnoty v týchto intervaloch budú zodpovedať farbám terra-in. colors(12). (1) 0 = (pn, //i2, ofi, af2, pi, p2i, p22, a2l, a 22, f>2)T ap = n\j(n\ +^2); parametre sú odhadnuté z dát. (2) 0 = (Jli, /í2, cr2, a2, p)T; parametre sú odhadnuté zo spoločného výberu. Vzorové riešenie pozri na obrázku 9 (dáta two-samples-correlations-trunk. txt). DU Odlišnosti od teoretického rozdelenia. Odlišnosti empirického rozdelenia (rozdelenia realizácií) od teoretického (napr. normálneho) rozdelenia, môžeme charakterizovať napr. ako pravostranné alebo ľavostranne zošikmené rozdelenie (obrázok 10, prvý riadok vľavo a vpravo), ploché alebo špicaté rozdelenie (obrázok 10, prvý riadok uprostred). Pri viacrozmerných rozdeleniach je situácia komplikovanejšia. Pri dvojrozmernom normálnom rozdelení može byť napr. zošikmená jedna alebo obe premenné (príklad zošikmenia oboch premenných zľava pozri na obrázku 10, dolný riadok). Príklad 23 (binomické rozdelenie, binomický experiment) Experiment pozostávajúci z fixného počtu Bernouliho experimentov (ozn. N) sa nazýva binomický experiment. Pravdepodobnosť úspechu ozn. p, pravdepodobnosť neúspechu q = 1 — p. Náhodná premenná X je počet pozorovaných úspechov počas experimentu. Pravdepodobnosť X = x za podmienky, že X pochádza z binomického rozdelenia Bin(N,p) píšeme ako Pr(X = x) = (^^(1 — p)N~x,x = 0,1,..., N (Ugarte a kol. 2008). Stredná (14. decembra 2014) Katina, S., 2014: Štatistická inferencia I 11 priemerná dlzka dolnej končatiny (mm) priemerná dlzka dolnej končatiny (mm) priemerná dlzka dolnej končatiny (mm) Obr. 9: Spoločná hustota dvojrozmerného normálneho rozdelenia (vľavo), hustota zmesi dvoch dvojrozmerných normálnych rozdelení (uprostred) a dvojrozmerný jadrový odhad superponovaný hustotou zmesi dvoch dvojrozmerných normálnych rozdelení (vpravo) - reálne dáta hodnota E [X] = N p a rozptyl Var[X] = Np(l—p). Naprogramujte a zobrazte v <Ä pravděpodobnostní! funkciu a (kumulatívnu) distribučnú funkciu pre Bin(5,0.5). Riešenie pozri na obrázku 11. cvič. Príklad 24 (multinomické rozdelenie) Majme náhodné premenné (1) socioekonomický status (vysoký - H, nízky - Lo), (2) politická príslušnosť (demokrat - D, republikán - R) a (3) politická filozofia (liberál - Li, konzervatívec - C). Označme ich interakcie nasledovne X\ (H-D-Li), X2 (H-D-C), X3 (H-R-Li), X4 (H-R-C), X5 (Lo-D-Li), X6 (Lo-D-C), X7 (Lo-R-Li) a X8 (Lo-R-C). Predpokladajme, že máme náhodný výber s rozsahom N = 50. Pravdepodobnosti p j sú nasledovné D-Li D-C R-Li R-C spolu H 0.12 0.12 0.04 0.12 0.4 Lo 0.18 0.18 0.06 0.18 0.6 spolu 0.30 0.30 0.10 0.30 1.0 Vypočítajte Var[X{\, Var[X/\, Cov [Xi,X/\, Cor [Xi,X/\ a očakávané početnosti N p j, j = 1,2,... ,8. pred Riešenie X = (XX,X2, ■ ■ ■ ,X8) ~ Mult(N,p), kde N = 50, p = (pi,p2,... ,Ps)T, vieme, že Xj ~ Bin(N,pj), p j sú v tabuľke v zadaní príkladu a j = 1,2,..., 8. Potom Var[Xľ] = 50 x 0.12 x (1 - 0.12) = 5.28, Var[Xs] = 50 x 0.04 x (1 - 0.04) = 1.92. Vybraná kovariancia a korelácia (medzi počtami príslušných skupín) je rovná Cov [Xi,X3] = -50 x 0.12 x 0.04 = -0.24, Cor [Xi,X3] = -0.24/V5.28 x 1.92 = -0.075. Očakávané počtetnosti pre každú bunku tabulky sú (všeobecne nemusia byť) celé čísla: D-Li D-C R-Li R-C H 6 6 2 6 Lo 9 9 3 9 (14. decembra 2014) Katina, S., 2014: Štatistická inferencia I 12 N(0,1) t{0,i,df=i0) sN(0,1) st{0,i,df=i0) -4 -2 0 2 4 - N{0,0.5) - N(0,1) - N{0,1.5) SN{0,0.5) - sN(0,1) - sN{0,1.5) I I -4 -2 I I I 0 2 4 Ll = 0, Oi = 0.5, o-i =1, o2 = 1.5 N(0,1) SN(0,1) -3-2-10 1 2 3 -3-2-10 1 2 3 4 -3-2-10 1 2 3 Ll-I = 0, |12 = 0, G-\ = 1, o2 = 1, p = 0.5 * "• •"' *jif8jäK' / "í \1^^ůnry* ^ 0.02 / I I I I I -3 -2-10 1 1 I 2 3 Ll-I = 0, |12 = 0, O-] = 1, 02 = 1 p = 0.5 Obr. 10: Hustoty normálneho rozdelenia a zošikmeného normálneho rozdelenia pri rôznych parametroch (prvý riadok); hustoty dvojrozmerného zošikmeného normálneho rozdelenia (druhý riadok vl'avo a uprostred) a dvojrozmerného normálneho rozdelenia (druhý riadok vpravo) pri rôznych parametroch Príklad 25 (súčinové multinomické rozdelenie) Majme dáta z predchádzajúceho príkladu a náhodný výber s rozsahom N± = 30 zo skupiny H, ďalší náhodný výber s rozsahom N2 = 20 zo skupiny Lo. Označme interakcie premenných nasledovne Xn = Xi\i (H-D-Li), X12 = X2\ľ (H-D-C), X13 = X3{1 (H-R-Li), X14 = Xtn (H-R-C), X21 = X1{2 (Lo-D-Li), X R-Li) a X24 - x = (x1;x2; s Xj\k, kde j -22 X2\2 (Lo-D-C), X23 = X3{2 (Lo--- X4{2 (Lo-R-C), kde X1 = (Xn,X12,X13,X14)T a X2 = (X21,X22,X23,X24)T. Potom má súčinové multinomické rozdelenie s K = 2, N± = 30, J\ = 4, N2 = 20,J2 = 4. Zápis = 1,2,3,4 a k = 1,2 zvýrazňuje fakt, že rozdelenie je podmienené socioekonomickým statusom (vysoký - H, nízky - Lo), t.j. rozdelenie v stĺpcoch tabulky je podmienené jej riadkom. Realizácie Xj\k označujeme ako rij\k = nkj, pravdepodobnosti ekvivalentné Xj\k = Xkj ako pj\k = pkj. Vypočítajte podmienené pravdepodobnosti pj\k, očakávané početnosti Nkpkj, Var[Xi3], Cov [X2i,X23] aCor[Xn,X23]. Riešenie Pravdepodobnosti štyroch kategórií asociovaných s H statusom sú podmienené pravdepodobnosti dané H statusom. Napr. Pr(X3|i) = 0.04/0.4 = 0.1. Pr(Xi|i) = 0.12/0.4 = 0.3, Pr(X3|2) = 0.06/0.6 = 0.1. Musíme ale tabulku prepísať na súčinovo-multinomický model, teda podmienené pravdepodobnosti pj\i dané socioekonomickým statusom i budú D-Li H 0.3 Lo 0.3 D-C 0.3 0.3 R-Li 0.1 0.1 R-C 0.3 0.3 spolu 1.0 1.0 pred (14. decembra 2014) Katina, S., 2014: Štatistická inferencia I 13 Q- ^ - pravděpodobnost™ funkcia pre X~Bin(5, 0.5) distribučná funkcia pre X~Bin(5, 0.5) Obr. 11: Pravdepodobnostná funkcia a distribučná funkcia Bin(5,0.5) Pre Ni = 30 a N2 = 20 máme očakávané počty nasledovné H Lo D-Li 9 6 D-C 9 6 R-Li 3 2 R-C spolu 9 30 6 20 Var(X3ll) = 30 x 0.1 x (1 - 0.1) = 2.7. Vybrané kovariancie (medzi počtami prislušných skupín) sú rovné Cov Cov Xl\l, X; 3|2 -20 x 0.3 x 0.1 = -0.6, 0, lebo Xi a X2 sú nezávislé. Príklad 26 (Poissonovo rozdelenie; počet havárií za týždeň) Ak každý z 50 miliónov ľudíšoféruje auto v Taliansku budúci týždeň nezávisle, potom pravdepodobnosť smrti pri autonehode bude 0.000002, kde počet úmrtí má binomické rozdelenie Bin(50mil, 0.000002) alebo limitné Poissonovo rozdelenie s parametrom 50mil x 0.000002 = 100. Príklad 27 (Poissonovo rozdelenie; pruské armádne jednotky) Nech početnosti úmrtí X ako následok kopnutia koňom v Pruských armádnych jednotkách (Bortkiewicz 1898), má Poissonovo rozdelenie s parametrom X, t.j. X ~ Poiss(X). Pravdepodobnosť, že niekto bude smrteľne zranený v danom dni je extrémne malá. Majme 10 vojenských jednotiek za 20-ročnú periódu s rozsahom M = 200 (200 = 10 x 20), kde popri početnostiach úmrtí n = 1,2,3,4,5+, v danej jednotke a v danom roku, zaznamenávame aj početnosti vojenských jednotiek mn pri danom n, kde M = £mn (pozri tabuľku). Vypočítajte očakávané početnosti, za predpokladu X ~ Poiss(X), kde X E„»m-y n 0 1 2 3 4 5+ mn 109 65 22 3 1 0 Príklad 28 (podiel chlapcov a dievčat v rodinách) Nech X predstavuje početnosť chlapcov medzi deťmi v rodinách. Tu môžeme predpokladať, že X ^ Bin(N,p), t.j. rodina môže mať vychýlený pomer pohlaví detí v smere ku chlapcom alebo dievčatám. V realite teda môžeme mať príliš veľa rodín len s chlapcami alebo len s dievčatami a nemáme dostatok rodín s pomerom pohlaví blízkym 51 : 49 (pomer chlapcov ku dievčatám). Z toho nám vyplýva, že rozptyl početnosti chlapcov bude v skutočnosti väčší ako rozptyl predpokladaný binomickým modelom Bin(N,p). cvic. pred (14. decembra 2014) Katina, S., 2014: Štatistická inferencia I 14 Príklad 29 (overdispersion v binomickom modeli) V klasickej štúdii pomeru pohlaví u ľudí z roku 1889 na základe záznamov z nemocníc v Sasku (bližšie pozri Lindsey a Altham (1998)) Geissler (1889) zaznamenal rozdelenie počtu chlapcov v rodinách. Medzi M = 6115 rodinami s N = 12 deťmi pozoroval nasledovné početnosti chlapcov (n sú početnosti chlapcov a mn početnosti rodín s n chlapcami) n 0 1 2 3 4 5 6 7 8 9 10 11 12 mn 3 24 104 286 670 1033 1343 1112 829 478 181 45 7 Vypočítajte mn za predpokladu, že početnosti chlapcov X v rodinách majú binomické rozdelenie s parametrami n = j^^nn = 0.5192 a N = 12, ozn. X ~ Bin(N ,-k). cvič. Riešenie n 0 1 2 3 4 5 6 7 8 9 10 11 12 očakávané mn 1 12 72 258 628 1085 1367 1266 854 410 133 26 2 Keď porovnáme pozorované mn a vypočítané (teoretické) mn zistíme, že pozorované poukazujú na overdispersion, t.j. máme väčšie početnosti rodín s malým a velkým množstvom chlapcov v porovnaní s teroretickými početnosťami. Príklad 30 (overdispersion v Poissonovom modeli) Majme početnosti úrazov n medzi robotníkmi v továrni, kde početnosti robotníkov mn pri danom n pozri v tabuľke (Greenwood a Yule 1920). cvič. n 0 1 2 3 4 > 5 mn 447 132 42 21 3 2 Vypočítajte očakávané početnosti robotníkov za predpokladu, že početnosti úrazov na robotníka X majú Poissonove rozdelenie s parametrom A = = 0.47, ozn. X ~ Poiss(X). Riešenie n 0 1 2 3 > 4 očakávané mn 406 189 44 7 1 Keď porovnáme pozorované mn a vypočítané (teoretické, očakávané) mn zistíme, že pozorované poukazujú na overdispersion, t.j. máme viac robotníkov bez úrazu ako aj viac robotníkov s väčším množstvom úrazov v porovnaní s teroretickými početnosťami. Príklad 31 (binomický rozdelenie, simulačná štúdia) Vygenerujte pseudonáhodné čísla X (početnosti úspechov) opakovane M-krát (M = 1000,) z Bin{N,p), kde N = 5 a p = 0.5. Vytvorte tabulku vygenerovaných (simulovaných) ako aj teoretických relatívnych početností (pre n = 0,1,... ,5). Superponujte histogram vygenerovaných pseudonáhodných čísel s teoretickou pravdepodobnostnou funkciou, cvič. Riešenie (pozri obrázok 12) r 0 1 2 3 4 5 simulované relatívne početnosti teoretické relatívne početnosti 0.021 0.031 0.152 0.156 0.318 0.312 0.324 0.312 0.155 0.156 0.030 0.031 (14. decembra 2014) Katina, S., 2014: Štatistická inferencia I 15 úspechy X Obr. 12: Histogram vygenerovaných pseudonáhodných čísel superponovaný spojnicovým grafom teoretickej pravdepodobnostnej funkcie X Príklad 32 (binomické vs normálne rozdelenie) Nech X n ~ Bin{N,p), potom môžeme aproximoval binomické rozdelenie normálnym nasledovne - Xn ~ N(Np, Np(l — p)), kde tiež platí Zn = ^y™~^P j ~ ^(0,1)- Ukážte, že CLV platí pre N = 100 a p = 1/2 na tri desatinné miesta. Riešenie Príklad hovorí o tom, ako dobre normálne rozdelenie aproximuje binomické pri rozsahu N = 100, čo je dôležité pri testovaní hypotéz. E[XN] = Np = 50, y/Var [XN] = y/Np(l-p) = y/E. Ak YN = XN/N, potom Pr^ - 1/2| < e) = 0.236, kde e = 0.02. Pr(0.48 < Y100 < 0.52) = Pr(48 < úpravy na spojitosť. Xim < 52) = Pr(48.5 < X100 < 51.5) = pr(^=^ < zim < ^50), kde Xim ~ N (50, 5) s použitím Príklad 33 (normálne rozdelenie, simulačná štúdia) Na základe simulačnej štúdie preverte, že ak X ~ iV(150,6.25), potom Xn ~ A^(150,6.25/n). Použite n = 30. Pre každú simuláciu X vypočítajte aritmetické priemery xm, m = 1,2,..., M, kde M = 500000. Superponujte ich histogram v relatívnej škále s teoretickou krivkou hustoty pre Xn. Vypočítajte Pr(Xn > 151) zo simulovaných dát a porovnajte tento výsledok s teoretickou (očakávanou) pravdepodobnosťou. Riešenie pozri na obrázku 13. cvič. Príklad 34 (normálne rozdelenie, simulačná štúdia) Nech X ~ N([ii,af) a Y ~ N(fi2,a2). _ _ 2 2 P otom X ril—Y ri2 ~ N(pi—fi2, ^+^)- Generujte pseudonáhodné čísla X aY rozdelení N'(// j, o~J), j = 1, 2, kde pi = 100, o~\ = 10, fi2 = 50, o2 = 9 pri (a) nľ = 4, n2 = 5, (b) nľ = 100, n2 = 81. Pre každú simuláciu X a Y vypočítajte rozdiel xm — ym,m = 1,2,..., M, kde M = 1000. Superponujte histogram týchto rozdielov v relatívnej škále s teoretickou krivkou hustoty rozdielu Xni — Yn2. Pre prípad (a) aj (b) vypočítajte Pr(Xni — Yri2) < 52 na základe empirického (vygenerovaného) a teoretického rozdelenia Xni — Yri2. cvič. Príklad 35 (štatistika) Majme náhodný výber (X\,X2,..., Xn)T, kde Xť 6 l,i = 1,2,...,«, potom príkladmi štatistík sú: Tx = Y/U Xí e r> T2 = Eľ=i Xi £R+U {0}, T3 = (YU Xh YU xľ )2 pred (14. decembra 2014) Katina, S., 2014: Štatistická inferencia I 16 Príklad 36 (testovacia štatistika, simulačná štúdia) Na základe simulačnej štúdie preverte, že ak náhodná premenná X má asymptoticky binomické rozdelenie Bin(N,p), potom testovacia štatistika Zw = x/n~p = > má~ asymptoticky normálne rozdelenie N(0,1). Použite p = 0,0.1,0.5,0.9 a 1, a V (p(ľ-p)/N N = 5,10,30,50 a 100. Okomentujte výsledky v spojitosti s Haldovou podmienkou Np(l — p) > 9. Pre každú simuláciu X vypočítajte z\y,m, m = 1,2,..., M, kde M = 1000. Superponujte histogram vygenerovaných testovacích štatistík v relatívnej škále s teoretickou krivkou hustoty Zy/. cvič. Príklad 36 hovorí o použití j ednovýberovej testovacej štatistiky pre parameter binomického rozdelenia (pravdepodobnosť) pre rôzne pravdepodobnosti a rôzne početnosti. Ak Haldova podmienka nie je splnená, nie je možné testovaciu štatistiku použiť. Príklad 37 (testovacia štatistika, simulačná štúdia) Na základe simulačnej štúdie preverte, že ak (a) X ~ N(fi,a2), kde = 0, a2 = 1 a (b) X ~ Exp(X), kde X = 1, E[X] = 1 a Var[X] = 1, potom testovacia štatistika F = ^n^S má asymptoticky Xn-i rozdelenie s n — 1 stupňami volnosti. Použite rozsahy náhodných výberov n = 15 a n = 100. Pre každú simuláciu X vypočítajte Fpoz>m, m = 1,2,..., M, kde M = 1000. Superponujte histogram vygenerovaných testovacích štatistík v relatívnej škále s teoretickou krivkou hustoty F. cvič. Príklad 38 (maximálne vierohodné odhady; Poissonovo rozdelenie) Každý rok za posledných päť rokov boli v nejakom meste registrované 3, 2, 5, 0 a 4 zemetrasenia za rok. Za predpokladu, že počet zemetrasení za rok X má Poissonovo rozdelenie s parametrom X, t.j. X ~ Poiss(X), odhadnite X (X predstavuje očakávanú početnosť zemetrasení za rok). pred Príklad 39 (T(p) a rozptyl pre p; X ~ Bin(N,p)) Z funkcie vierohodnosti odvoďte pozorovanú Fisherovu mieru informácie X(p) a rozptyl Var\p\. pred Príklad 40 (X(X) a rozptyl pre A; X ~ Poiss(X)) Každý rok za posledných päť rokov boli v nejakom meste registrované 3, 2, 5, 0 a 4 zemetrasenia za rok. Za predpokladu, že počet zemetrasení za rok X ~ Poiss(X), odhadnite rozptyl parametra X a vypočítajte hodnotu tohoto odhadu rozptylu pre počet zemetrasení. pred Príklad 41 (kvadratická aproximácia funkcie vierohodnosti) (1) Nakreslite škálovaný logaritmus funkcie vierohodnosti binomického rozdelenia. Na x-ovej osi bude p a na y-ovej osi ln£(p) = í(p|x) — max(/(p|x)). Porovnajte ln£(p) s kvadratickou aproximáciou vypočítanou pomocou Taylo-rovho rozvoja \n£(p) = ln(^|ľ|^) « — |X(p)(p — p)2. (2) Nech skóre funkcia S (p) = JMnL(p|x). Keď zoberieme deriváciu kvadratickej aproximácie uvedenej vyššie, dostaneme S (p) ~ — T (p) (p — p) alebo —T~1l2(p)S(p) ~ X1/2(p)(p — p). Potom zobrazením pravej strany na x-ovej osi a ťavej strany na y-ovej osi dostaneme asymptoticky lineárnu funkciu s jednotkovým sklonom. Asymptoticky tiež platí X1^2(p)(p — p) ~ N(0,1). Je postačujúce mať rozsah x-vej osi (—2,2), pretože funkcia je asymptoticky (lokálne) lineárna na tomto intervale. Rozumne škálujte y-ovú os. Zobrazte pre (a) n = 8, N = 10, (b) n = 80, N = 100 a (c) n = 800, N = 1000 (p E (0.5,0.99);. Okomentujte rozdiely medzi (a), (b) a (c). Grafické riešenie je na obrázku 13. cvič. Príklad 42 (1(0) pre vektor 0 = (/x,o-2)T; X ~ N(fi, a2)) Nech X ~ N(fi,a2). Čomu je rovná pozorovaná Fisherova informačná matica X(0), kde 0 = (fl,32)T ? pred (14. decembra 2014) Katina, S., 2014: Štatistická inferencia I 17 1 i i i i i i i i i i i i i r -2-1012 -2-1012 -2-1012 skalovana 9 skalovana 9 skalovana 9 (a) linearita skoré funkcie (b) linearita skoré funkcie (c) linearita skoré funkcie Obr. 13: Porovnanie škálovaného logaritmu funkcie vierohodnosti (plná čiara) s jeho kvadratickou aproximáciou (čiarkovaná čiara) v prvom riadku a porovnanie škálovanej skóre funkcie a priamky s nulovým interceptom a jednotkovým sklonom v druhom riadku Príklad 43 (X(p) a rozptyl pre p; X ~ Multj(N,p)) Z funkcie vierohodností odvoďte pozorovanú Fisherovu informačnú maticu X(p) a kovariančnú maticu Var\p\. pred Príklad 44 (profilová vierohodnosť; normálne rozdelenie) Profilová funkcia vierohodnosti pre H počítaná pre každé fixované fi, kde maximálne vierohodný odhad a2 bude a2 = ^ Yľi=i(xi ~ A4)2' má tvar L(/x|x) = c (a2 Yn/2, kde c je nejaká konštanta. L(fi\x.) nie je identická s odhadnutou funkciou vierohodnosti L(fi,a2 = íx2|x) = cexp (—^ Ylr=i(xi ~ f1)2)> t-J- s rezom L(fi,a2\x.) v bode a2 = d2. Obe funkcie vierohodností budú veľmi podobné, ak je rozptyl a2 dobre odhadnutý. V opačnom prípade sa preferuje profilová funkcia vierohodností. Profilová funkcia vierohodností pre a2 je rovná L(a2\x) = c (a2yn/2 exp (-^ £ľ=ifc - x)2) = c (a2)~n/2 exp (-na2/(2a2)). pred Príklad 45 (kvadratická aproximácia profilovej funkcie vierohodnosti) (1) Nakreslite šká-lovaný logaritmus profilovej funkcie vierohodnosti normálneho rozdelenia pre fi. Na x-ovej osí bude u a na y-ovej osí ln£(/x|x) = Z(//|x) — max(/(p|x)). Porovnajte ln£(/x|x) s kvadratickou aproximáciou vypočítanou pomocou Taylorovho rozvoja ln£(/x|x) = ln(^|^|^) « — |X(/x)(/x — /i)2. (2) Nech skóre funkcia S(jjl) = ^MnL(//|x). Keď zoberieme deriváciu kvadratickej aproximácie uvedenej vyššie, dostaneme S(/i) ~ — X(/x)(/x — //) alebo —X^1/2('J1)S(jjl) X1/2(/í)(/í — //). Potom zobrazením pravej strany na x-ovej osi a ľavej strany na y-ovej osi dostaneme asymptoticky lineárnu funkciu s jednotkovým sklonom. Asymptoticky tiež platí X1/2(/i) (p — //) ~ N(0,1). Je postačujúce mai rozsah x-vej (14. decembra 2014) Katina, S., 2014: Štatistická inferencia I 18 osi (—2,2), pretože funkcia je asymptoticky (lokálne) lineárna na tomto intervale. Rozumne školujte y-ovú os. Zobrazte pre (a) n = 10, (b) n = 100 a (c) n = 1000. Použite (í) X ~ N(0,1) a (2) X ~ (1 — p)N(0,1) + pN(0,2), kde p = 0.05. Okomentujte rozdiely medzi (a), (b) a (c), ako aj rozdiely medzi (í) a (2). DU Príklad 46 (maximálne vierohodný odhad fi a a2) Vygenerujte pseudonáhodné čísla z X ~ N (A, 1), n = 1000. (a) Napíšte logaritmus profilovej funkcie vierohodnosti pre fi a a2 a preverte, či sú maximálne vierohodné odhady fi a a2 dostatočne blízko k ich skutočným hodnotám. Nakreslite grafy í(/x|x) a /(o"2|x); kde zvýrazníte polohu maxím týchto funkcií, (b) Napíšte logaritmus funkcie vierohodnosti pre 0 = (fi, a2) a preverte, či je maximálne vierohodný odhad 0 = (fi,a2)T dostatočne blízko k jeho skutočnej hodnote, (c) Nakreslite graf /((//, o-2) |x) použitím funkcie image() a superponujte ho s konturovým grafom použitím funkcie contourO. Zvýraznite polohu maxima (pozri obrázok 14). DU V o2 V Obr. 14: Profilová funkcia vierohodnosti pre fi (vľavo), a2 (uprostred) a funkcia vierohodnosti pre oba parametre (vpravo); X ~ N (A, 1)); maximálne vierohodné odhady strednej hodnoty a rozptylu sú označené zvislou čiarkovanou čiarou (vľavo a uprotred) a maximálne vierohodný odhad vektora parametrov je označený • (vpravo) Príklad 47 (binomické rozdelenie, maximálne vierohodný odhad p) Nech X ~ Bin(N,p) a realizácie X sú x = n. Predpokladajme, že sme pozorovali (a) x = 2, (b) x = 10 a (c) x = 18 úspechov v N = 20 pokusoch. Pomocou <® vypočítajte maximálne vierohodný odhad p. Výsledok zobrazte do grafu spolu s funkciou vierohodnosti. cvič. Príklad 48 (maximálne vierohodné odhady; multinomické rozdelenie) Majme dáta more-samples-probabilities-pubis. txt. Nakreslite logaritmus štandardizovanej funkcie vierohodnosti v parametroch pľ ap2 Európskej populácie (nľ = 30, n2 = 20 a n3 = 10) pomocou funkcie contourO. Dokreslite do obrázku jej maximum v bode 0 = (pi,P2)T- DÚ Príklad 49 (overdispersion v Poissonovom modeli, pokrač.) Majme početnosti úrazov n medzi robotníkmi v továrni, kde početnosti robotníkov mn pri danom n pozri v tabuíke (Greenwood a Yule 1920). n 0 i 2 3 4 > 5 mn 447 132 42 21 3 2 Vypočítajte mn za predpokladu, že početnosti úrazov na robotníka X majú negatívne binomické rozdelenie s parametrami a a -k. cvič. (14. decembra 2014) Katina, S., 2014: Štatistická inferencia I 19 2 Charakteristiky polohy a variability a štatistická grafika Príklad 50 (argument minima) Vygenerujte pseudonáhodné čísla X ~ N(fi,a2), n = 1000, fi = 0, a2 = 1. Vygenerované čísla ozn. xi} i = 1,2,..., 1000. Nájdite numericky také c, ktoré minimalizuje (a) sumu štvorcov odchýlok ^^°(^i — c)2, t.j. c\ = arg mirivc Y^í^=l (xí ~ c)2 a (b) sumu absolútnych odchýlok \xí — c\, t.j. c2 = arg minvc Y^í^=l \xí — c\. Za c dosadzujte postupne (1) všetky x^ (x^ sú usporiadané Xi podia velkosti od najmenšieho po najväčšie) a vybrané charakteristiky polohy ako (2) aritmetický priemer, (3) nejaké kvantily xp, kde p G (0,1) a pod. Nakreslite obrázok závislosti (a) sumy štvorcov odchýlok na x^, t.j. body [x j, y j], kde y j = ^^°(^i — x^)2 a (b) sumy absolútnych odchýlok na x^, t. j. body [x^,yj], kde y j = Y^l^i \xí ~ x(j)\- Podobné obrázky nakreslite aj pre xp namiesto x^y DU Definícia 1 (asymptotické rozdelenie poriadkovej štatistiky) Nech X^, X(2), ■ ■ ■, X(n) sú poriadkové štatistiky náhodného výberu X\, X2,..., Xn. Majme pravdepodobnosť a, kde F(ta) = a. Asymptoticky platí, že \fň(^ — a) konverguje k 0. Potom je poriadková štatistika X^ normálne rozdelená so strednou hodnotou E[X^] = ta a rozptylom crx = "2^""^■ Ak X ~ N(fi,a2), potom °X{i) = ° 241nn 1 Príklad 51 (rozptyl poriadkovej štatistiky) Pomocou delta metódy odvoďte rozptyl poriadkovej štatistiky v definícii 1. DU Príklad 52 (rozptyl poriadkovej štatistiky, X ~ N(fi,a2)) Pomocou definície 1 odvoďte rozptyl poriadkovej štatistiky, ak X ~ N(fi,a2). DU Definícia 2 (stredná hodnota a rozptyl mediánu) Stredná hodnota mediánu X^n+i^ je rovná E[X(n+i^] = Jí a rozptyl mediánu o\ = 4p^n, kde n je nepárne. Ak X ~ N(fi,a2), potom ^(H+i) =°2^- Pred Príklad 53 (rozptyl mediánu) Pomocou delta metódy odvoďte rozptyl poriadkovej štatistiky v definícii 2. DÚ Príklad 54 (rozptyl mediánu, X ~ N (u, a2)) Pomocou definície 2 odvoďte rozptyl poriadkovej štatistiky, ak X ~ N(fi,a2). DÚ Príklad 55 (pásy normality) Na základe vygenerovaných pseudonáhodných čísel X ~ N(fi,a2), n = 100, fi = 0,a2 = 1, kde M = 1000, odhadnite (a) hustotu m-tej realizácie pomocou funkcie densityO; ponechajte argument n=512 a nastavte from=-3 a to=3; (b) distribučnú funkciu m-tej realizácie pomocou (a) a funkcie cumsumO a (c) empirické kvantily m-tej realizácie pomocou funkcie qqnorm(). Vygenerované čísla xmi,m = 1,2,..., 1000 a i = 1,2, ...,100, uložte po riadkoch do matice X, ktorá bude mať rozmery 1000 x 100. Odhadnuté hustoty a distribučné funkcie uložte po riakodch do matíc H a D, ktoré bude mať rozmery 1000 x 512 a empirické kvantily do matice K, ktorá bude mať rozmery 1000 x 100. Pre každú z matíc H, D a K vypočítajte 2?o.o5 a xo.95 P° stĺpcoch a zobrazte ich ako pásy pomocou funkcie polygonO. Do obrázkov vkreslíte (a) teoretickú hustotu, (b) teoretickú distribučnú funkciu a (c) kvantilovú priamku (pomocou funkcie qqlineO) červenou farbou. Obrázky usporiadajte ako trojicu vedľa seba. Dáta, ktorých normalitu chceme graficky testovať budú (1) X ~ N{0,1), n = 100, (2) X ~ [pN(0,1) + (1 - p)N(0,4)], n = 100 a p = 0.95 a (3) X ~ [pN(0,1) + (1 - p)iV(0,4)], n = 100 a p = 0.9. Zobrazte (í), (2) a (3) oddelene do grafov (a), (b) a (c). Okomentujte. Riešenie (t) pozri na obrázku 15. (14. decembra 2014) Katina, S., 2014: Štatistická inferencia I 20 t-1-1-1-1-1-r-^ "n-1-1-1-1-1-r-^ ^-1-1-1-1-1-r -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3 data data teoretické kvantily Obr. 15: Pásy spoľahlivosti normálneho rozdelenia - pre hustotu (vľavo), distribučnú funkciu (uprostred) a kvantilovú priamku (vpravo) 3 Testovanie hypotéz Aby sme mohli z dát vyvodiť nejaké interpretovateľné závery, musíme najprv ciele (účely, biologicky formulované hypotézy) výskumu preformulovať do matematicko-štatistickej podoby. Takáto formulácia je dôležitá pri výbere správneho modelu na dáta. V prípade parametrického modelu máme jeho parametre 9, kde hypotézy sú o týchto parametroch. Jednou hypotézou je tzv. nulová hypotéza, čo je tvrdenie o 9, ktoré testujeme prostredníctvom nejakého štatistického testu. Ďalšou je alternatívna hypotéza, ktorá je doplnkom nulovej hypotézy v priestore, z ktorého pochádza 9. Na základe výsledku testu vyvodíme záver o 9. Proces testovania hypotéz nazývame aj štatistická inferencia, ktorej základom je štatistická teória testovania hypotéz. Definícia 3 (štatistický test) Ak parametrický priestor O, do ktorého patri parameter 9 modeluj-' rozdelíme na dve podmnožiny ©o a ©i, kde O = ©oU©i, potom štatistický test T = T(X) predstavuje funkciu X z výberového priestoru V do priestoru 0O U ©i. Statistický test je teda pravidlo výberu jednej z dvoch možnosti, nulovej a alternatívnej, na základe dát. Nulovú hypotézu definujeme ako Hq : 9 G ©o a alternatívnu hypotézu ako H\ : 9 G ©i. Prvky množiny V, pre ktoré zamietame Hq, predstavujú tzv. oblasť (obor) nezamietania nulovej hypotézy. Tieto prvky označujeme 3^o ■ V tomto prípade Hq nezamietame. Prvky množiny V, pre ktoré nezamietame Hq, predstavujú tzv. oblasť (obor) zamietania nulovej hypotézy alebo kritickú oblasť. Tieto prvky označujeme y± = Wx- V tomto prípade Hq zamietame. Rozdelenie priestoru 0 na obor zamietania WT = W a nezamietania At prebieha na základe testovacej štatistiky T (X). Stretávame sa s nasledovnými štyrmi možnosťami: A) ak platí Hq, tak naše rozhodutie je nezamietnuť Hq (správne): B) ak platí Hq, tak naše rozhodutie je zamietnuť Hq (nesprávne): C) ak neplatí Hq, tak naše rozhodutie je nezamietnuť Hq (nesprávne): D) ak neplatí Hq, tak naše rozhodutie je zamietnuť Hq (správne). V prípade B, C je naše rozhodnutie nesprávne, v prípade A, D správne. V prípade B sa dopúšťame chyby prvého druhu (CHPD), kde Pr(CHPD) < a a a sa nazýva hladina významnosti. Doplnok (14. decembra 2014) Katina, S., 2014: Štatistická inferencia I 21 ku a je pravdepodobnosť 1 — a a nazýva sa koeficient spoľahlivosti (prípad A). V prípade C sa dopúšťame chyby druhého druhu (CHDD), kde Pr(CHDD) = f3. Doplnok ku Pr(CHDD) je pravdepodobnosť 1 — f3 a nazýva sa sila testu (pri nejakej alternatíve; prípad D). Sila testu závisí na zvolenej testovacej metóde a hlavne na tom, aké je skutočné rozdelenie náhodnej premennej X. aký je typ použitej testovacej štatistiky alebo aké sú skutočné hodnoty parametrov 9. Vyššie uvedené pravdepodobnosti môžeme zosumarizovať nasledovne: A) 1 — a < Pr(nezamietni Hq\Hq platí) = Pr(nezamietni Hq\Hi neplatí): B) a > Pr(CHPD) = Pr(zamietni Hq\H0 platí) = Pr(zamietni Hq\Hi neplatí): C) f3 = Pr(CHDD) = Pr(nezamietni Hq\H0 neplatí) = Pr(nezamietni Hq\Hi platí): D) 1 — f3 = Pr(zamietni Hq\Hq neplatí) = Pr(zamietni Hq\Hi platí). V prípade (A) a (B) hovoríme o teste na hladine významnosti a (a-level test). Ak v (A) a (B) nahradíme znamienko nerovnosti rovnosťou, hovoríme o teste úrovne a (size a test). Test T (x) charakterizujeme jeho silofunkciou f3*(9) = 1 — f3 (9) = Pľ(9 G 0 : T(x) = ©i). Táto je závislá na 9 a niekoľkých iných argumentoch (vždy sú špecifikované pri konkrétnych testoch) ako aj na type alternatívnej hypotézy. Hladinu významnosti a je možné definovať tiež pomocou f3*(9), kde a = sup8eQ0 P*(9), t.j. ide o najväčšiu možnú pravdepodobnosť chyby I. druhu. Hladina významnosti a je daná (určená štatistikom, experimentátorom) vopred, testovacia štatistika T(x) sa vypočíta na základe realizácií (odpovedí) x náhodnej premennej X, kritická hodnota (alebo kvantil) sa nájde v štatistických tabulkách alebo ju vypočítame pomocou *®. V praxi teda voľba W závisí na požiadavke, aby pravdepodobnosť chyby I. druhu bola menšia alebo rovná zvolenému kladnému číslu a, kde a G (0,1/2), najčastejšie a = 0.05, 0.01 alebo 0.001. Súčasne volíme W tak, aby pravdepodobnosť chyby I. druhu bola čo najmenšia. Optimálne by sme chceli mať f3*(9) (pri nejakej hodnote 9) tak velké číslo (blížiace sa jednotke) ako je možné, keď 9 G ©i a také malé číslo, keď 9 G ©o- Tieto dve požiadavky sú v konflikte, t.j. • zväčšovanie oboru nezamietania Hq zmenšuje pravdepodobnosť chyby I. druhu, ale zväčšuje pravdepodobnosť chyby II. druhu: • zmenšovanie oboru nezamietania Hq zväčšuje pravdepodobnosť chyby I. druhu, ale zmenšuje pravdepodobnosť chyby II. druhu. Extrémnymi prípadmi sú dve situácie, kde Tq(x) = ©o a Tí (x) = ©i, t.j. každý z týchto dvoch testov je najlepší možný, keď 9 patrí určitej podmnožině 0, ale najhorší možný v opačnom prípade. Klasickým prístupom k tejto dileme je Neyman-Pearsonov prístup, kde dopredu fixujeme a a vyberieme taký test, ktorý má pri danej alternatíve najväčšiu silu 1 — f3 pre 9 G ©i medzi všetkými možnými na teoreticky stanovenej (danej, nominálnej) hladine významnosti a. Testy spomínané v kapitolách 5 až 8 sú najsilnejšími testami pri daných alternatívach za určitých predpokladov (bližšie pozri kap. 5 až 8). V reálnych situáciách (alebo často aj za porušenie predpokladov testov) sa nominálna hladina významnosti líši od aktuálnej (skutočnej) hladiny významnosti. Teda rozhodnutie / skutočnosť Hq nezamietnuť Hq zamietnuť Hq platí správne rozhodnutie chyba I. druhu Hq neplatí chyba II. druhu správne rozhodnutie (14. decembra 2014) Katina, S., 2014: Štatistická inferencia I 22 Definícia 4 (kvantil) Nech Fx je distribučná funkcia náhodnej premennej X a a G (0,1), kde Fx {xa) = a. Potom číslo xa = F^1 (a) sa nazýva a-kvantil príslušného rozdelenia. Pre vybrané rozdelenia platí • Pr (X < x(l — a)) = 1 — a pre X ~ N (0,1) (normálne rozdelenie s parametrami fi = 0 a a2 = 1; hovorí sa mu aj normálne normované rozdelenie alebo štandardizované normálne rozdelenie); x(l — a) sa často zapisuje ako x\_a; • Pr [X < Xdf(l ~ a)) = 1 — ct pre X ~ x% (chi-kvadrát rozdelenie s df stupňami volnosti); • Pr (X < tdf (1 — a)) = 1 — a pre X ~ tdf (Studentovo t-rozdelenie s df stupňami voľnosti); • Vr (X < Fdf14f2{^ — ex)) = 1 — a pre X ~ Fdf1,df2 (Fisherovo F-rozdelenie s df\ a df2 stupňami voľnosti). Hustoty vyššie spomenutých rozdelení pozri na obrázku 16. Tiež platí • Pr (xa/2 < X < ari_a/2) = Fx (^i-a/2) - Fx (xa/2) = 1 — a, kde a E (0,1/2): • Pr (X < Qi) = 1/4, kde Qx = Fxľ (1/4) sa nazýva výberový prvý kvartil (dolný kvartil): . Pr (X < Q2) = 1/2, kde Q2 = Fxľ (1/2) sa nazýva výberový druhý kvartil (medián): • Pr (X < Q3) = 3/4, kde Q3 = F^1 (3/4) sa nazýva výberový tretí kvartil (horný kvartil): • Qd = F^1 {kj 10) sa nazýva výberový /c-ty decil, Qp = F^1 {kj 100) sa nazýva výberový fc-ty percentil. d - N(0,1) ....... t(5) .......... t(100) chi-sq(5) chi-sq(10) chi-sq(20) -4 -2 0 10 20 30 40 50 Obr. 16: Grafické znázornenie hustôt normálneho rozdelenia, í-rozdelenia, %2-rozdelenia a F-rozdelenia pri rôznych stupňoch voľnosti Definícia 5 (kritická hodnota) Kritická hodnota príslušného rozdelenia je hodnota, ktorú náhodná premenná X prekročí s pravdepodobnosťou a. Pre vybrané rozdelenia platí • Pr (X > u(a)) = a pre X ~ N (0,1); u(a) sa často zapisuje ako ua; (14. decembra 2014) Katina, S., 2014: Štatistická inferencia I 23 • Pr (X > x%{ tdf (a)) = a pre X ~ tdf; • Pr (X > Fdfudf2 (a)) = a pre X ~ Fdfudh. Z vyššie uvedeného je zreteľné, že a x 100% kritická hodnota je identická s (1 — a) x 100% kvantilom. Pri kritickej hodnote počítame obsah pod krivkou hustoty príslušného rozdelenia nad týmto bodom a pri kvantile pod týmto bodom. Stupne voľnosti df predstavujú množstvo nezávislej informácie, ktoré potrebujeme na charakterizáciu rozdelenia pravdepodobnosti nejakej náhodnej premennej. Používajú sa aj namiesto n ako prevencia pred nadhodnotením odhadu nejakého parametra (napr. rozptylu). Ide najčastejšie o rozsah náhodného výberu n zmenšený o jedna (napr. Studentovo í-rozdelenie, F-rozdelenie), kde jednotka odpočítaná od n znamená jeden voľne viazaný parameter (napr. jedna stredná hodnota, jeden rozptyl a pod.). Pri zložitejších modeloch s viacerými parametrami odpočívame od n počet voľne viazaných parametrov (napr. počet stredných hodnôt) alebo df môže tiež predstavovať počet hladín kategoriálnej premennej zmenšený o počet voľne viazaných parametrov. Príklad 56 (štandardizované normálne rozdelenie) Vypočítajte kritické hodnoty u(a) rozdelenia N(0,1), kde a = 0.1, 0.05, 0.01, 0.025 a 0.005. Riešenie v X(n). Majme transformáciu T{1) = Fn(X{1)), T(2) = Fn(X(2)),...,T(n) = Fn(X{n)). Potom T{1), T(2), ..., T(n) stí poriadkové štatistiky. Potom platí lim Pr(sup [Fn(ar) - F(x)]n1/2 < A) = $(A), n->oo V:re-y fcde -F(X) je teoretická distribučná funkcia a $(A) = J^fcl_00(—l)fce~2fc2A2. Potom 100 x (1 — «)% pas spoľahlivosti pre Fn(x) definujeme ako Fn(x) ± A^l/n1/2, Ä;de $(Aa) = 1 — a a Var[Fn(x)] = l/n. Potom môžeme tvrdiť, že F (X) patri do 100 x (1 — a)% pásu spoľahlivosti a zároveň je medzi nulou a jednotkou s pravdepodobnosťou 1 — a. pred Príklad 62 (graf distribučnej funkcie a jej IS) Nakreslite graf distribučnej funkcie N(fi,a2), kde fi = 0 a a2 = 1. Do grafu dokreslite 95% pás spoľahlivosti pre F (x). Jeho hranice vypočítajte pomocou simulácie pseudonáhodných čísel z N(0,1) pri n = 50, kde Fn(x) je odhadnutá z dát. cvič. (14. decembra 2014)