Maja Pohar Perme Verjetnost in statistika z nalogami Ljubljana, 2014 Skripte Ekonomske fakultete Maja Pohar Perme Verjetnost in statistika z nalogami Založila : Ekonomska fakulteta v Ljubljani, Založništvo Šifra: POH14AU114 Recenzenta: dr. Simona Korenjak Černe dr. Marko Pahor Objavljeno na spletni strani: http://www.ef.uni-lj.si/zaloznistvo/studijska_gradiva CIP – Kataložni zapis o publikaciji Narodna univerzitetna knjižnica, Ljubljana 519.2(0.034.2) POHAR Perme, Maja Verjetnost in statistika z nalogami [Elektronski vir] / Maja Pohar Perme. - El. knjiga. - Ljubljana : Ekonomska fakulteta, 2014 ISBN 978-961-240-275-4 (pdf) 272503296 Vse pravice pridržane. Noben del gradiva se ne sme reproducirati ali kopirati v kakršni koli obliki: grafično, elektronsko ali mehanično, kar vključuje (ne da bi bilo omejeno na) fotokopiranje, snemanje, skeniranje, tipkanje ali katere koli druge oblike reproduciranja vsebine brez pisnega dovoljenja avtorja ali druge pravne ali fizične osebe, na katero bi avtor prenesel materialne avtorske pravice. Kazalo 1 Verjetnost 1 1.1 Normalna porazdelitev . . . . . . . . . . . . . . . . . . . . . . 2 1.2 Generiranje spremenljivk . . . . . . . . . . . . . . . . . . . . . 8 1.3 Vsota diskretnih slučajnih spremenljivk . . . . . . . . . . . . . 12 1.4 Vsota zveznih slučajnih spremenljivk . . . . . . . . . . . . . . 14 1.5 Vsota normalnih spremenljivk . . . . . . . . . . . . . . . . . . 18 1.6 Porazdelitev vzorčnega povprečja . . . . . . . . . . . . . . . . 22 1.7 Pogojna pričakovana vrednost in varianca . . . . . . . . . . . 24 1.8 Uvrščanje . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 1.9 Centralni limitni izrek . . . . . . . . . . . . . . . . . . . . . . 35 1.10 Normalna aproksimacija . . . . . . . . . . . . . . . . . . . . . 37 2 Vzorčenje 41 2.1 Vzorčenje - neskončna populacija . . . . . . . . . . . . . . . . 42 2.2 Vzorčenje - končna populacija . . . . . . . . . . . . . . . . . . 45 2.3 Določitev načrta vzorčenja . . . . . . . . . . . . . . . . . . . . 48 2.4 Ocena kovariance . . . . . . . . . . . . . . . . . . . . . . . . . 51 2.5 Enostavni vzorec, končna populacija II . . . . . . . . . . . . . 55 2.6 Večstopenjsko vzorčenje . . . . . . . . . . . . . . . . . . . . . 58 3 Ocenjevanje parametrov - metoda največjega verjetja 64 3.1 Ocenjevanje deleža . . . . . . . . . . . . . . . . . . . . . . . . 65 3.2 Povezanost dveh spremenljivk . . . . . . . . . . . . . . . . . . 69 4 Preizkušanje domnev 74 4.1 Preizkušanje domnev . . . . . . . . . . . . . . . . . . . . . . . 75 4.2 Moč testa . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 4.3 Enostavni domnevi . . . . . . . . . . . . . . . . . . . . . . . . 81 i Verjetnost in statistika Kazalo 4.4 Enostavni domnevi, posplošitev . . . . . . . . . . . . . . . . . 87 4.5 Posplošeni test razmerja verjetij . . . . . . . . . . . . . . . . . 88 5 Linearna regresija 95 5.1 Linearna regresija . . . . . . . . . . . . . . . . . . . . . . . . . 95 5.2 Matrično računanje . . . . . . . . . . . . . . . . . . . . . . . . 99 5.3 Predpostavke linearne regresije . . . . . . . . . . . . . . . . . 105 ii Uvod Statistika je veda z zelo raznolikimi področji uporabe, vsako področje za- znamuje svoja narava podatkov in z njimi povezani zapleti. Na prvi pogled statistika lahko deluje le kot nepregleden skupek formul in receptov, ključ do njenega razumevanja se skriva v razumevanju skupnih temeljev statističnega sklepanja, ki so osnova vsem metodam in pristopom. Knjig na temo statističnega sklepanja je veliko, nekatere se skozi osnovne ideje statistike uspejo prebiti skorajda brez formul, druge so precej bolj teoretično zahtevne in za razumevanje potrebujejo veliko matematičnega znanja. Knjiga, ki jo imate v rokah, izbira srednjo pot, za razumevanje ni potrebno znanje teorije mere, še več, večino primerov je moč rešiti z rahlo obogate- nim znanjem srednješolske matematike. Navkljub temu je izpeljav veliko in marsikateremu nematematiku se bo snov vsaj na začetku zdela le težko ra- zumljiva. V nasprotju z večino knjig, ki so namenjene podajanju, izpeljevanju in doka- zovanju teorije in so nujen spremljevalec tega gradiva, želi ta knjiga študenta voditi skozi osnovne zakonitosti statistike s pomočjo podrobno rešenih primerov. Knjiga torej ni mišljena kot samostojno gradivo, ki podaja vse definicije in izreke, temveč bolj kot spremljevalec, ki študenta vodi po osnovah teoretične statistike in mu pomaga razumeti posamezne korake izpeljav. Pred- vsem pa snovi skuša dajati poudarke in študenta opominjati na to, zakaj so na videz zelo teoretične lastnosti tako pomembne za praktično uporabo statistike. Jedro knjige predstavlja 26 nalog, ki s svojimi rešitvami študentu odkrivajo osnovne zakonitosti in ideje statistike. Pomemben delež nalog najprej za- jame nujne osnove verjetnosti, nato sledijo povsem statistične teme. Vsako poglavje se prične s kratkim uvodom, ki pomaga umestiti naloge v snov, vsem rešitvam so dodani povzetki, ki študenta še enkrat opozorijo na bistvene po- udarke posamezne naloge. Primeri so namenoma izbrani z različnih področij iii in želijo študenta ves čas opominjati tako na široko uporabnost osnovnih sta- tističnih znanj kot tudi na možnost, da za nek problem, na katerega naletijo v okviru podatkov s področja financ, prav lahko že obstaja rešitev na nekem povsem drugem področju, na primer v biostatistiki. Skoraj vse naloge v besedilu so povsem originalne, a ker so zaradi jasnosti izpeljav in razumevanja izbrani dokaj osnovni primeri, je seveda marsikatero zelo podobno nalogo moč srečati tudi v drugih statističnih knjigah. Po- dane rešitve nalog niso edine možne in pravilne, študent je vedno vabljen, da najprej poišče svojo rešitev, nato pa s pomočjo danih rešitev razume tudi alternativne možnosti, predvsem pa glavne poudarke, ki iz rešitve sledijo. Za bolj samostojno reševanje pred vpogledom v rešitve so nekaterim nalogam dodani tudi namigi. Knjiga je namenjena tako matematikom kot tudi tistim z mnogo manj ma- tematičnega znanja. Za slednje prebijanje iz vrstice v vrstico pogosto pov- zroča precej frustracij, zato je na mestu opozorilo, da naj se ne obremenjujejo preveč s tehničnimi podrobnostmi, temveč skušajo raje razumeti splošne cilje in ugotovitve posamezne naloge. To opozorilo naj pride prav tudi tistim, ki se jim teorija zdi sorazmerno lahka - za pravo razumevanje statistike ne za- dostuje le podrobno razumevanje izpeljav, nujno je razumeti, kaj neki dokaz oziroma lastnost pomeni za rešitev nekega realnega problema. Knjiga po organizaciji snovi sledi knjigi Mathematical Statistics and Data Analysis (Rice, 2009), ki je zaradi podobne matematične zahtevnosti in so- sledja snovi vsekakor priporočen vir teorije pri reševanju nalog. Za lažjo orientacijo so pri posameznih snoveh navedena tudi poglavja te knjige, kjer je predstavljena potrebna teorija. Prav tako knjiga vsebuje veliko referenc na dodatno literaturo, s pomočjo katere študenti lahko poglobijo posamezne snovi. Zelo koristna vira sta na primer tudi Casella and Berger (1990) in Shao (2003), vendar je predvsem slednji precej bolj matematično zahteven. Pomembna lastnost statistične teorije, ki izhaja iz realnih problemov, je, da si lahko vse zakonitosti ogledamo tudi v praksi, z računalniškimi simulacijami. Za resnično razumevanje pojmov, kot so porazdelitev, varianca, nepristran- skost, standardna napaka, ipd., je zato neprecenljiva tudi zmožnost, da si jih študent prikaže s pomočjo ustreznega računalniškega programa. Skoraj vsem nalogam so na koncu zato dodani še predlogi, kako teorijo spoznavati s konkretnimi primeri v ustreznem statističnem programu. Za primere je izbran R (R Development Core Team, 2013), statistični program, ki predsta- vlja eno najbolj uporabljanih orodij med statistiki. Seveda je vse tovrstne naloge mogoče reševati tudi s poljubnim drugim programskim orodjem, ki iv omogoča generiranje podatkov. Rešitve teh nalog so prepuščene študentu, dodani so le kosi kode, ki naj bi študentu prihranili nekaj časa. Tako je pri prvih nalogah dodanih nekaj zelo osnovnih nasvetov in komentarjev kode, pri kasnejših nalogah pa je koda podana le, kadar je nekoliko zahtevnejša. Poglavje 1 Verjetnost Eden osnovnih ciljev statistične analize je, da na podlagi podatkov, zbra- nih na nekem, ponavadi sorazmerno majhnem vzorcu, skušamo sklepati o splošnih zakonitostih, ki veljajo v celotni populaciji. Zaradi naključne variabilnosti je na podlagi vzorca o splošnih zakonitostih seveda nemogoče karkoli trditi z gotovostjo, statistika zaupanje v sklepe o populaciji ovrednoti z verjetnostjo. Verjetnost je torej osnovno orodje statističnega sklepanja in prvo poglavje je namenjeno spoznavanju osnovnih pojmov in postopkov, ki jih bomo potrebovali v kasnejših, bolj ‘statističnih’ delih knjige. Če želimo govoriti o verjetnosti posameznih dogodkov, moramo poznati porazdelitev slučajne spremenljivke, ki opisuje možne izide. To poglavje se ukvarja z več porazdelitvami, ki se v statistiki pogosto pojavljajo (glej npr. Rice, 2009, razdelka 2.1 in 2.2), in raziskuje njihove osnovne lastnosti. Jedro poglavja predstavljajo štiri teme: izpeljava gostote transformirane spremenljivke (Rice, 2009, razdelek 2.3), formula za vsoto slučajnih spremenljivk (Rice, 2009, poglavje 3), pričakovana vrednost in varianca (Rice, 2009, po- glavje 4) ter centralni limitni izrek (Rice, 2009, poglavje 5). Prva naloga raziskuje lastnosti normalne porazdelitve, ki je prav gotovo najpomembnejši primer porazdelitve v statistični teoriji, ter pokaže, da s poljubno linearno transformacijo normalne spremenljivke spet dobimo normalno porazdelitev. Sledi naloga, ki predstavi enakomerno porazdelitev in z njo povezan pomem- ben statistični izrek (‘transformacija integrala verjetnosti’), ki je ključen pri generiranju slučajnih spremenljivk s poljubno drugo porazdelitvijo. V tretji nalogi izpeljemo formulo za vsoto dveh diskretnih slučajnih spremen- ljivk ter jo skušamo razumeti na preprostem primeru. Nato v četrti nalogi 1 Verjetnost in statistika 1.1. Normalna porazdelitev pokažemo, da lahko na povsem enak način razumemo tudi formulo za zve- zni primer. V peti nalogi formulo za vsoto uporabimo na primeru normalne spremenljivke in pokažemo, da tudi s seštevanjem neodvisnih normalnih spre- menljivk ostajamo pri normalni porazdelitvi. Poglavje se v šesti in sedmi nalogi nadaljuje s spoznavanjem lastnosti priča- kovane vrednosti, variance in kovariance, pri tem je izpostavljeno predvsem razumevanje formul za njihove pogojne vrednosti. Uporabnost teh spoznanj je nato prikazana pri izpeljavi lastnosti na videz povsem smiselnega načina uvrščanja, za katerega se izkaže, da deluje drugače, kot bi si želeli. Poglavje se zaključi s centralnim limitnim izrekom, ki bo igral bistveno vlogo v zakulisju marsikatere metode v nadaljevanju - za kakršnokoli statistično sklepanje bomo morali poznati porazdelitev, kadar ne bo znana oziroma bo njena izpeljava prezahtevna, jo bomo s pomočjo centralnega limitnega izreka skušali aproksimirati z normalno. 1.1 Normalna porazdelitev Krvni doping je metoda, pri kateri si športnik kri najprej odvzame, nato pa si jo včrpa pred pomembnim nastopom in tako umetno poveča število rdečih krvničk ter si s tem izboljša trenutno počutje in vzdržljivost. Ker ni udeleženih tujih substanc, krvnega dopinga ni mogoče neposredno odkriti. Zato ga skušajo odkrivati s statističnimi metodami - doping naj bi naka- zovale vrednosti krvnih parametrov (hemoglobina), ki pretirano narastejo (včrpanje) oz. padejo (odvzem). Vemo, da je vrednost hemoglobina pri nedopingiranem športniku porazde- ljena normalno s povprečjem µ = 148 in varianco σ2 = 85. Označimo vred- nost hemoglobina z X, torej X ∼ N (148,85). a) Izračunajte verjetnost, da je posameznikova vrednost večja od 166. V ta namen izpeljite formulo: • Naj bo X ∼ N (µ, σ2), kako je porazdeljena porazdeljena slučajna spremenljivka Y = aX + b, kjer je a > 0? Namig: Zapišite najprej porazdelitveno funkcijo, nato izrazite go- stoto. Ali lahko gostoto zapišete kot gostoto normalne spremen- ljivke? 2 Verjetnost in statistika 1.1. Normalna porazdelitev FY (y) = P (Y ≤ y) = P (aX + b ≤ y) = P (aX ≤ y − b) y − b y − b = P (X ≤ ) = FX( ) a a 1 y − b fY (y) = fX( ) a a Za normalno porazdeljeno X torej velja: 1 [ y−b − µ]2 f a Y (y) = √ exp − a · 2πσ2 2σ2 1 [y − (b + aµ)]2 = exp − 2π(a · σ)2 2(a · σ)2 Torej, Y ∼ N (a · µ + b, (a · σ)2). Linearna transformacija normalne spremenljivke je še vedno normalna. • Kaj moramo vzeti kot a in b, da bo Y standardno normalno poraz- deljena spremenljivka, torej Y ∼ N (0,1)? a mora biti enak 1 , b pa − µ . Uporabiti moramo torej transforma- σ σ cijo Y = X−µ in nato verjetnosti odčitati iz tabel za standardno σ normalno porazdelitev (oz. uporabiti ustrezno numerično metodo). V našem primeru je X = 166 in zato Y = X−148 √ = 166−148 √ = 1,95. 85 85 Iz tabel za standardno normalno porazdelitev (ali pa s pomočjo računalnika) izvemo, da je P (X ≤ 166) = P (Y ≤ 1,95) = 0,974, zato je verjetnost P (X > 166) = 0,026. b) Izračunajte (simetrične) meje, ki jih nedopingiran športnik preseže z ver- jetnostjo manj kot 0,01. Naj bo Y standardno normalno porazdeljena spremenljivka, zanimajo nas meje, izven katerih je vrednost te spremenljivke z verjetnostjo 0,01. Če želimo postaviti simetrične meje, pomeni, da nas zanimata tisti vrednosti, izven katerih je v repih na vsaki strani verjetnost 0,005. Iz tabel izvemo, da je P (Y ≥ 2,58) = 0,995, ustrezna mejna vrednost standardno nor- malno porazdeljene spremenljivke je torej ±2,58. 3 Verjetnost in statistika 1.1. Normalna porazdelitev Y = X−148 √ , zato 85 X − 148 X − 148 0,99 = P ( √ ≤ 2,58) + P ( √ > −2,58) 85 85 √ √ = P (X ≤ 148 + 2,58 · 85) + P (X > 148 − 2,58 · 85) = P (X ≤ 171,8) + P (X > 124,2) c) Naj bodo meje take, kot ste jih izračunali v prejšnji točki. Športnika te- stiramo 10x na leto. Kolikšna je verjetnost, da vsaj enkrat preseže meje (pri tem predpostavimo, da so meritve narejene v dovolj velikih časovnih presledkih, da so med seboj neodvisne)? Naj bo U Bernoullijevo porazdeljena spremenljivka U ∼ Ber(0,01), kjer je {U = 1} = {vrednost je izven meja}. Imamo 10 neodvisnih realizacij te slučajne spremenljivke, Ui, i = 1, . . . ,10, za vsako velja P (Ui = 1) = 0,01. Ker so neodvisne, velja P (U1 = 0, U2 = 0, . . . , U10 = 0) = {P (U1 = 0)}10. Verjetnost, da v 10 meritvah ne preseže meja, je 0,9910, verjetnost, da jih preseže vsaj enkrat, je P = 1 − 0,9910 = 0,096. d) Naj bo Y ∼ N (0,1). Izračunajte porazdelitev slučajne spremenljivke Y 2. Katero znano porazdelitev dobite? Namig: Porazdelitev gama ima gostoto fT (t) = λαtα−1e−λt . Γ(α) √ √ FZ(z) = P (Z ≤ z) = P (Y 2 ≤ z) = P (− z ≤ Y ≤ z) √ √ = FY ( z) − FY (− z) 1 √ 1 √ fZ(z) = √ fY ( z) + √ fY (− z) 2 z 2 z 1 √ √ = √ [fY ( z) + fY (− z)] 2 z 1 1 = √ [e−z/2 + e−z/2] = √ e−z/2 2 z · 2π z · 2π Gama porazdelitev ima gostoto fT (t) = λαtα−1e−λt . Če vzamemo, da je Γ(α) √ α = 1 in λ = 1 , ter upoštevamo, da je Γ( 1 ) = π, dobimo natanko gornjo 2 2 2 formulo. Torej je Y 2 ∼ Γ( 1 , 1 ) (to je hkrati tudi porazdelitev χ2). 2 2 1 4 Verjetnost in statistika 1.1. Normalna porazdelitev Opombi: • Pri zveznih porazdelitvah je verjetnost, da spremenljivka zavzame točno določeno vrednost P (Y = y), neskončno majhna, zato je vse- eno, ali v porazdelitveni funkciji uporabljamo znak ≤ ali <, torej P (Y ≤ y) = P (Y < y). • Kadar želimo pokazati, da lahko gostoto spremenljivke zapišemo kot gostoto neke že znane porazdelitve, je dovolj pokazati enako funk- cijsko obliko argumentov. Konstanta bo potem vedno enaka, saj je integral gostote po celem definicijskem območju vedno enak 1. e) Raziskovalci na področju športa so dokazali, da je pri biatloncih hemo- globin izven tekmovalnega obdobja porazdeljen kot N (150, 80), med tek- movalnim obdobjem pa kot N (146, 80). Tekmovalno obdobje je pri teh športnikih dolgo približno pol leta. Zanima nas porazdelitev hemoglo- bina, če ne vemo, kdaj je bil vzorec odvzet. Ali je ta porazdelitev še vedno normalna? Namig: P (A) = n P (A|B B i=1 i)P (Bi), kadar velja P ( n i=1 i) = 0 in P ( n B i=1 i) = 1. Definiramo Bernoullijevo porazdeljeno spremenljivko Y , ki naj označuje obdobje (0 = izven, 1 = tekme, verjetnost vsakega izida je 0,5). Poznamo pogojni porazdelitvi: Z|Y = 0 ∼ N (150, 80), Z|Y = 1 ∼ N (146, 80). Porazdelitev Z je torej (uporabimo namig, kjer je B1 = {Y = 0} in B2 = {Y = 1}, namesto z verjetnostmi pišemo z gostotami) fZ(z) = fZ|Y =0(z)P (Y = 0) + fZ|Y =1(z)P (Y = 1) 1 1 = fZ|Y =0(z) + fZ|Y =1(z) 2 2 1 = √ e− (z−146)2 2·80 [1 + e− 16−8(z−146) 2·80 ] 2 2π80 Ta spremenljivka v splošnem ni normalno porazdeljena, glej sliko 1.1. Predlogi za vaje v R • Generirajte 10000 realizacij normalno porazdeljene spremenljivke X ∼ N (148,85) (rnorm). Narišite histogram (hist), izračunajte delež vrednosti nad 166 (sum(x>166)/10000). 5 Verjetnost in statistika 1.1. Normalna porazdelitev 400 300 enca 200 Frekv 100 0 100 120 140 160 180 Z Slika 1.1: Porazdelitev spremenljivke Z (za ta primer je za povprečje v tek- movalnem obdobju vzeta vrednost 120, da bi bila porazdelitev bolj različna od normalne). • Oglejte si funkcijo qnorm in z njo poiščite meje, izven katerih je športnik z verjetnostjo 0,01. Primerjajte z deležem v vašem primeru. • Transformirajte vrednosti spremenljivke X tako, da dobite standardno normalno porazdeljeno spremenljivko (y=(x-148)/sqrt(85)). Preverite grafično s histogramom. • Generirajte po 10 vrednosti za 10000 posameznikov. Izračunajte delež posameznikov, ki imajo vsaj eno vrednost izven intervala [124,2, 171,8]. • Narišite histogram za vrednosti X2, primerjajte z rezultatoma funkcij rgamma in rchisq. • Narišite še porazdelitev slučajne spremenljivke iz zadnje točke (upora- bite bolj različni povprečji, da se prepričate, da porazdelitev zares ni normalna). > set.seed(1) #dolocimo seme, da bodo rezultati vedno enaki > z0 <- rnorm(1000,mean=150,sd=sqrt(80)) #1000 vredn.izven tek.obdobja > z1 <- rnorm(1000,mean=120,sd=sqrt(80)) #1000 vredn. v tek.obdobju > z <- c(z0,z1) #obe obdobji enako pogosto zastopani v novi spr. > hist(z,main="",xlab="Z",ylab="Frekvenca") #prikaz porazdelitve 6 Verjetnost in statistika 1.1. Normalna porazdelitev Povzetek • Pomembna lastnost, ki precej olajša računanje z normalno porazdeli- tvijo, je, da je linearna transformacija normalne spremenljivke zopet normalno porazdeljena. To je tudi razlog, zakaj je dovolj, da so stati- stične tabele vedno podane le za ‘standardno normalno porazdelitev´ - z njeno pomočjo lahko namreč preprosto izračunamo verjetnosti za normalno porazdeljeno spremenljivko s poljubnima vrednostima para- metrov. Seveda pa navkljub tej lepi lastnosti ne smemo pričakovati, da bo kar vsaka transformacija normalne porazdelitve ostala normalna - protiprimer je podan v točki (e) te naloge. • V nalogi je prikazan postopek izpeljevanja porazdelitve spremenljivke Y , ki je definirana kot neka transformacija spremenljivke X z znano po- razdelitvijo. Kumulativno porazdelitveno funkcijo nove spremenljivke FY najprej izrazimo z znano kumulativno porazdelitveno funkcijo FX, nato pa uporabimo, da je gostota fY odvod funkcije FY . • Točka (c) je dodana kot opozorilo, da so vse izračunane verjetnosti ve- ljavne le pri enem preizkusu. Čim bomo nek preizkus uporabili večkrat, ne smemo pozabiti na problem večkratnega preizkušanja , torej dejstva, da se verjetnosti pri večjem številu opravljenih preizkusov močno spre- menijo. • Točka (e) je v to nalogo umeščena predvsem kot vaja kako formalno zapisati nek besedni opis porazdelitve. Dani zapis pri tem seveda ni edini možen, je pa iz njega jasno vidno, da ne gre za normalno poraz- delitev. V tej točki smo uporabili pogojno porazdelitev, ki je v knjigi Rice (2009) predstavljena v razdelku 3.5. • V nalogi smo izpeljali še rezultat, da je kvadrat standardno normalno porazdeljene spremenljivke spremenljivka porazdeljena z gama poraz- delitvijo: 1 1 X ∼ N (0,1) ⇒ X2 ∼ Γ( , ). 2 2 Gama porazdelitev s tema parametroma je hkrati enaka tudi porazde- litvi χ2 z eno stopinjo prostosti. 7 Verjetnost in statistika 1.2. Generiranje spremenljivk 1.2 Generiranje slučajnih spremenljivk s po- močjo enakomerne porazdelitve Generator (psevdo)slučajnih vrednosti iz enakomerne spremenljivke zgene- rira željeno število vrednosti xi, ki so porazdeljene kot X ∼ U [0, 1]. a) Kako bi s pomočjo tega generatorja dobili 10 realizacij Bernoullijevo po- razdeljene spremenljivke Y , pri kateri je P (Y = 1) = 0,1? Generirajmo 10 vrednosti, npr.: > set.seed(4) > runif(10) [1] 0.585800305 0.008945796 0.293739612 0.277374958 [5] 0.813574215 0.260427771 0.724405893 0.906092151 [9] 0.949040221 0.073144469 Vrednostim, ki so pod 0,1, damo vrednost 1, ostalim pa 0, torej: > set.seed(4) > (runif(10)<0.1)*1 [1] 0 1 0 0 0 0 0 0 0 1 b) Recimo, da imamo spet 10 enot, vendar jim želimo dati različne verje- tnosti, da bodo izžrebane. Prvih pet enot želimo izžrebati z verjetnostjo 0,3, drugih pet pa z verjetnostjo 0,1 (kot primer si zamislimo žreb, v ka- terem želimo dati prednost ženskam. Verjetnost za vsakega posameznika v našem vzorcu določimo glede na spol - prvih pet je žensk, drugih pet je moških). Kako bi iz istim generatorjem zagotovili ustrezno porazdelitev? > set.seed(4) > (runif(10) 0). Kako bi jih lahko simulirali z uporabo prej omenjenega genera- torja? Najprej potrebujemo funkcijo F : z FZ(z) = λe−λxdx 0 λ z = e−λx −λ 0 = −1[e−λz − 1] = 1 − e−λz Inverzna porazdelitev F −1 je enaka: u = 1 − e−λx 1 − u = e−λx − log(1 − u) = λx − log(1 − u) x = λ Če so vrednosti u torej realizacije enakomerno porazdeljene slučajne spre- menljivke U , so x realizacije eksponentno porazdeljene spremenljivke X. Mimogrede - log tu in v celotnem preostalem besedilu označuje naravni logaritem, torej logaritem z osnovo e. f) Kako bi hkrati simulirali vrednosti za posameznike z različno vrednostjo λ? Enako kot zgoraj - le da so vrednosti λ lahko različne. Generiranje eksponentno porazdeljenih vrednosti bo osnova za generiranje časov pri analizi preživetja - osnovna predpostavka na tem področju so namreč enakomerno porazdeljeni časi od diagnoze neke bolezni do smrti. Parameter eksponentne porazdelitve je pri tem lahko drugačen za vsakega posameznika, odvisen je namreč lahko npr. od njegove starosti, spola, ipd. (glej zadnji primer v R). 10 Verjetnost in statistika 1.2. Generiranje spremenljivk Predlogi za vaje v R • Generirajte podatke za 10000 voznikov, tako da jih je 500 med njimi pijanih, 9500 pa treznih. Naj bo verjetnost, da ima pijani voznik av- tomobilsko nesrečo 0,3, verjetnost za zdravega pa 0,003. Izračunajte delež nesreč na simuliranih podatkih in ga primerjajte z dejansko ver- jetnostjo nesreče. • Generirajte podatke za 100 posameznikov, tako da bo njihova starost enakomerno porazdeljena med 50 in 80. Preverite s histogramom. • Vzemimo, da je bila posameznikom iz prejšnje točke postavljena dia- gnoza hude bolezni. Generirajte čase preživetja z eksponentno poraz- delitvijo, tako da bodo imeli starejši posamezniki večjo verjetnost, da umrejo prej. Namig: parameter λ naj bo premosorazmeren s starostjo, npr. λ = starost/100. Povzetek • Pomemben izrek (transformacija integrala verjetnosti ), ki ga izpeljemo v tej nalogi, pove, da če spremenljivko X preslikamo z njej lastno ku- mulativno porazdelitveno funkcijo FX, vedno dobimo enakomerno po- razdeljeno spremenljivko. Poleg izpeljave dokaza je dodan še poskus intuitivne razlage. Pokažemo, da velja tudi obratno - če vrednosti neke enakomerno porazdeljene spremenljivke preslikamo z inverzom neke po- razdelitvene funkcije F −1, bo F Y Y kumulativna porazdelitvena funkcija novo dobljene spremenljivke. Tako lahko s pomočjo generatorja ena- komerne porazdeljenih slučajnih spremenljivk dobimo poljubno drugo porazdelitev, za katero lahko zapišemo inverz F −1 (ta seveda ne obstaja vedno - protiprimer najdemo že pri normalni porazdelitvi, kjer ekspli- citna formula za integral gostote in s tem kumulativno porazdelitveno funkcijo ne obstaja). • Izrek je v nalogi podan za dva primera porazdelitve - Bernoullijevo in eksponentno. Generiranje poljubne porazdelitve bo prišlo zelo prav v nalogah, ki sledijo - le če znamo generirati podatke iz znanih porazdeli- tev, lahko teoretično izpeljane lastnosti preverimo tudi s simulacijami. 11 Verjetnost in statistika 1.3. Vsota diskretnih slučajnih spremenljivk 1.3 Vsota diskretnih slučajnih spremenljivk Naj bosta X in Y neodvisni Bernoullijevo porazdeljeni spremenljivki, X, Y ∼ Ber(p). a) Kako je porazdeljena njuna vsota? Označimo Z = X + Y . Verjetnost, da je P (Z = z) za nek z, zapišemo kot vsoto verjetnosti vseh kombinacij X = x in Y = y, ki dajo vsoto enako z (lahko seštevamo, saj so dogodki nezdružljivi): P (Z = z) = P (X + Y = z) = P (X = z − y, Y = y) y = P (X = z − y|Y = y)P (Y = y) y Za neodvisni X in Y torej velja P (Z = z) = P (X = z − y, Y = y) = P (X = z − y)P (Y = y) y y Uporabimo, da sta spremenljivki Bernoullijevo porazdeljeni, in dobimo: P (Z = 0) = P (X + Y = 0) = P (X = 0)P (Y = 0) = (1 − p)2 P (Z = 1) = P (X = 0)P (Y = 1) + P (X = 1)P (Y = 0) = (1 − p)p + p(1 − p) = 2p(1 − p) P (Z = 2) = P (X = 1)P (Y = 1) = p2 Velja torej 2 P (Z = z) = pz(1 − p)2−z z b) Kako je porazdeljena vsota n neodvisnih enako porazdeljenih Bernoulli- jevih spremenljivk? Označimo Z = n X i=1 i, i = 1, . . . , n, Xi ∼ Ber(p). Za n = 2 je Z porazdeljena binomsko, Z ∼ Bin(2, p). Dokazati želimo, da to velja za 12 Verjetnost in statistika 1.3. Vsota diskretnih slučajnih spremenljivk poljuben n. V ta namen bomo uporabili matematično indukcijo - prvi ko- rak smo že pokazali, sledi še korak n → n + 1. Naj bo torej U ∼ Bin(n, p) in X ∼ Ber(p), pokazati moramo, da velja U + X ∼ Bin(n + 1, p). Za 1 ≤ z ≤ n velja P (Z = z) = = P (U = z − x)P (X = x) x = P (U = z)P (X = 0) + P (U = z − 1)P (X = 1) n n = pz(1 − p)n−z · (1 − p) + pz−1(1 − p)n−z+1 · p z z − 1 n n = + pz(1 − p)n−z+1 z z − 1 n + 1 = pz(1 − p)n+1−z z Dokaz za z = 0 in z = n + 1 je prepuščen bralcu. Predlogi za vaje v R • S pomočjo funkcije sample generirajte po 100 realizacij dveh Bernoul- lijevih spremenljivk in si oglejte porazdelitev njune vsote. Povzetek • Naloga služi kot uvod v naslednjo nalogo - na najpreprostejšem pri- meru diskretnih spremenljivk smo intuitivno skušali razumeti formulo za porazdelitev vsote dveh slučajnih spremenljivk: vsota dveh spremenljivk je lahko enaka z na več različnih načinov, in ker slučajna spremenljivka ne more zavzeti dveh različnih vrednosti hkrati, lahko dogodek X + Y = z razbijemo na posamezne nezdružljive dogodke (za vsako vrednost Y = y posebej). Njihove verjetnosti lahko zaradi nezdružljivosti enostavno seštejemo. Verjetnost vsote dveh diskretnih slučajnih spremenljivk lahko računamo s formulo P (X + Y = z) = P (X = z − y, Y = y) y 13 Verjetnost in statistika 1.4. Vsota zveznih slučajnih spremenljivk Če sta spremenljivki neodvisni, se formula poenostavi v P (X + Y = z) = P (X = z − y)P (Y = y) y • V nalogi smo izpeljali še rezultat, da je vsota n neodvisnih in enako porazdeljenih Bernoullijevih spremenljivk porazdeljena z binomsko po- razdelitvijo: n Xi ∼ Ber(p), i = 1 . . . n, Xi neodvisne ⇒ Xi ∼ Bin(n, p) i=1 1.4 Vsota zveznih slučajnih spremenljivk Poleg posamičnih vrednosti želijo pri športnikih proučevati tudi zaporedje več meritev. Zanima nas porazdelitev vsote kvadriranih standardiziranih odmikov od povprečja pri ničelni domnevi, da športnik ni kriv. Naj bo torej Z standardizirani odmik od povprečja (po predpostavki normalno porazdeljen), zanima nas Z2 (gledamo vsoto kvadriranih odmikov, saj so vrednosti lahko negativne ali pozitivne). Pri tem predpostavimo, da so bile meritve narejene v dovolj dolgih časovnih presledkih, da so vrednosti med seboj neodvisne. a) Najprej nas zanima, kako je porazdeljena vsota dveh neodvisnih zve- znih spremenljivk (izpeljite formulo za dve zvezni spremenljivki, torej Z = X + Y , primerjajte jo s formulo za diskretne). Zapišemo ustrezno kumulativno porazdelitveno funkcijo kot integral več- razsežne porazdelitve v ustreznih mejah: ∞ z−y P (Z ≤ z) = P (X + Y ≤ z) = fX,Y (x, y)dxdy −∞ −∞ ∞ z = fX,Y (v − y, y)dvdy −∞ −∞ pri čemer smo naredili substitucijo x = v − y. Sedaj lahko integrala zamenjamo in odvajamo (privzamemo, da je zunanji integral zvezen v z) 14 Verjetnost in statistika 1.4. Vsota zveznih slučajnih spremenljivk ter dobimo: z ∞ FZ(z) = fX,Y (v − y, y)dydv −∞ −∞ ∞ fZ(z) = fX,Y (z − y, y)dy −∞ ∞ fZ(z) = fX(z − y)fY (y)dy −∞ Pri tem zadnja vrstica velja, če sta X in Y neodvisni slučajni spremen- ljivki. Dobili smo rezultat, ki je analogen diskretni verziji. b) Kako se porazdeljuje S = Z2 + Z2, če sta spremenljivki Z 1 2 1 in Z2 poraz- deljeni standardno normalno in med seboj neodvisni? Namig: Pri izpeljavi uporabite rezultat a 1 dx = π 0 (a − x)x Izpeljali smo že, da je Z2 ∼ χ2 oziroma Z2 ∼ Γ( 1 , 1 ), torej 1 2 2 1 z f √ Z2 (z) = exp − ; z > 0 2πz 2 Izračunajmo gostoto vsote Z2 + Z2. Za dano vrednost s vemo, da mora 1 2 biti vrednost z2 med 0 (z1 in z2 ne moreta biti negativni) in s (ker sta obe pozitivni in je njuna vsota enaka s), zato morajo biti meje integracije 15 Verjetnost in statistika 1.4. Vsota zveznih slučajnih spremenljivk med 0 in s. s fS(s) = fZ2(s − z2)fZ2(z2)dz2 1 2 0 1 s 1 s − z 1 z = √ exp − 2 √ exp − 2 dz2 2π 2 z 2 0 s − z2 2 1 s s 1 1 = exp − √ √ dz2 2π 2 z 0 s − z2 2 1 s = exp − π 2π 2 1 s = exp − , 2 2 Za izračun integrala smo pri tem uporabili v namigu podano formulo. Gostota gama porazdelitve je λaxa−1e−λx fX(x) = ; x > 0, λ > 0, a > 0, Γ(a) Dobljeni rezultat je torej porazdelitev gama z λ = 1 in a = 1. 2 Na enak način bi lahko izpeljali tudi splošnejši rezultat: X1 ∼ Γ(a1, λ), X2 ∼ Γ(a2, λ) ⇒ X1 + X2 ∼ Γ(a1 + a2, λ) c) Denimo, da so športnikovi standardizirani odmiki (vrednosti Z) na petih merjenjih naslednji: 1,6; 1,5; −1,6; 1,8; 1,4. Kaj lahko sklepamo? Uporabimo, da je vsota n neodvisnih enako porazdeljenih spremenljivk Xi ∼ Γ( 1 , 1 ) porazdeljena kot n X , 1 ). 2 2 i=1 i ∼ Γ( n 2 2 > vr <- c(1.6,1.5,-1.6, 1.8, 1.4) > rez <- sum(vr^2) > rez [1] 12.57 > 1-pgamma(rez, 2.5, 0.5) [1] 0.02775943 16 Verjetnost in statistika 1.4. Vsota zveznih slučajnih spremenljivk Naša ničelna domneva je, da športnik ni kriv. Pod to ničelno domnevo se vsota kvadriranih odmikov porazdeljuje po gama porazdelitvi Γ( 5 , 1 ). 2 2 Verjetnost, da je vsota 12,57 ali več, je približno 0,03. Predlogi za vaje v R • Generirajte 10 vrednosti iz porazdelitve X ∼ N (148,85), te vrednosti naj predstavljajo 10 meritev pri enem športniku. Vrednosti standardi- zirajte (tako da dobite spremenljivko porazdljeno kot N (0,1)), kvadri- rajte in seštejte. To naj bo vrednost za prvega športnika, na enak način generirajte vrednosti za 1000 športnikov. Narišite histogram vredno- sti. S pomočjo funkcije pgamma poiščite mejo, ki jo porazdelitev Γ( 10 , 1 ) 2 2 preseže z verjetnostjo manj kot 0,01, in izračunajte delež športnikov na vašem vzorcu, ki presežejo to mejo. • Recimo, da ima dopingiran športnik enako povprečje, a večjo vari- anco (vrednosti bolj nihajo, saj manipulira s krvjo). Generirajte 1000 športnikov z večjo varianco in si oglejte, kolikšen delež bo presegel meje iz prejšnje točke. Povzetek • Formula za vsoto dveh zveznih slučajnih spremenljivk je analogna for- muli za vsoto diskretnih spremenljivk - verjetnosti zamenjajo gostote, vsote pa integrali: gostoto vsote dveh zveznih slučajnih spremenljivk lahko zapišemo kot ∞ fX+Y (z) = fX,Y (z − y, y)dy −∞ oziroma če sta spremenljivki neodvisni, kot ∞ fX+Y (z) = fX(z − y)fY (y)dy −∞ • V nalogi smo izpeljali še rezultat, da je vsota dveh neodvisnih gama porazdeljenih slučajnih spremenljivk (z enakim parametrom λ) zopet 17 Verjetnost in statistika 1.5. Vsota normalnih spremenljivk gama porazdeljena slučajna spremenljivka. Podobno lastnost - da je vsota neodvisnih enako porazdeljenih spre- menljivk zopet lahko opisana z isto porazdelitvijo - srečamo tudi pri normalni porazdelitvi (to dokažemo v nalogi 1.5), seveda pa ta lastnost še zdaleč ne velja za poljubno verjetnostno porazdelitev, na primer za vsoto Bernoullijevo ali enakomerno porazdeljenih spremenljivk (glej tudi povzetek naloge 1.9). 1.5 Vsota normalnih spremenljivk Pokazali smo že, da je linearna transformacija normalno porazdeljene spre- menljivke zopet normalno porazdeljena. V tej vaji bomo pokazali, da je vsota dveh neodvisnih (standardno) normalno porazdeljenih spremenljivk zo- pet normalno porazdeljena. Nato bomo s protiprimerom pokazali, da to ne velja za vsoto poljubnih odvisnih normalno porazdeljenih spremenljivk. a) Naj bosta spremenljivki X ∼ N (0,1) in Y ∼ N (0,1) in med seboj neod- visni. Kako je porazdeljena njuna vsota? Uporabimo formulo za gostoto vsote dveh neodvisnih slučajnih spremen- ljivk: ∞ fZ(z) = fX(z − y)fY (y)dy −∞ ∞ 1 (z − y)2 1 y2 = √ exp − √ exp − dy 2π 2 2π 2 −∞ ∞ 1 z2 = exp − exp −y2 exp {zy} dy 2π 2 −∞ Sedaj dele izraza, v katerih nastopa y, zapišemo kot kvadrat neke vsote: z z 2 z 2 y2 − zy = y2 − 2y + − 2 2 2 z 2 z2 = y − − 2 4 18 Verjetnost in statistika 1.5. Vsota normalnih spremenljivk Gornji integral torej prepišemo v ∞ 1 z2 z 2 z2 fZ(z) = exp − exp − y − exp dy 2π 2 2 4 −∞ ∞ 1 z2 (y − z )2 = exp − exp − 2 dy 2π 4 2 · 12 −∞  ∞  1 z2 1 1 (y − z )2 = √ exp − √ exp − 2 dy   2π 4 2 2π 1 2 · 12 −∞ 2 1 z2 = √ exp − 2π · 2 2 · 2 V predzadnji vrstici smo pod integralom dobili ravno gostoto normalne porazdelitve N ( z , 1 ), njen integral je zato 1 (ne glede na vrednost z, ki je 2 2 znotraj tega integrala konstanta). Spremenljivka Z je normalno porazde- ljena, Z ∼ N (0,2). b) Naj bosta X in Y neodvisni standardno normalno porazdeljeni spremen- ljivki, Z pa enaka |Y |, če je X ≥ 0, in Z = −|Y |, če je X < 0. Kako je porazdeljena spremenljivka Z? Najprej narišimo simulirane vrednosti z R: > set.seed(1) > x <- rnorm(1000,0,1) #1000 realizacij normalne spr., > y <- rnorm(1000,0,1) #povprecje=0, sd=1 > z <- abs(y) #z = |y| > z[x<0] <- -z[x<0] #z =-|y|, ce je x<0 > hist(z,main="",ylab="Frekvenca") #histogram z Sedaj še izpeljimo porazdelitveno funkcijo. Naj bo z < 0 (torej X nega- tiven), uporabimo, da sta spremenljivki X in Y neodvisni: FZ(z) = P (Z ≤ z) = P (X < 0, −|Y | ≤ z) = P (X < 0)[P (Y ≤ z) + P (Y ≥ −z)] 1 = [2 · P (Y ≤ z)] 2 = P (Y ≤ z) = FY (z) 19 Verjetnost in statistika 1.5. Vsota normalnih spremenljivk 150 enca 100 Frekv 50 0 −2 0 2 4 z Slika 1.3: Porazdelitev spremenljivke Z. Na enak način izpeljemo še P (0 ≤ Z ≤ z) = P (0 ≤ Y ≤ y) za z > 0. Pokazali smo, da je porazdelitvena funkcija Z enaka porazdelitveni funkciji Y , Z je torej standardno normalno porazdeljena spremenljivka. c) Skicirajte skupno porazdelitev spremenljivk X in Z. Ali sta spremenljivki neodvisni? Očitno je (glej sliko 1.4), da spremenljivki nista neodvisni, vedno imata enak predznak. ● ● ● ● 3 ● ● ● 3 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● 2 ● ● ● ● ● ● ● ● ● 2 ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ●● ● ● ● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● 1 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ●● ● ● ● ●● ● ● ● 1 ● ● ● ● ●● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●● ● ● ●● ● ● ● ● ●● ● ● ●● ●● ●● ● ● ●● ● ● ● ● ● ● ●● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● Y ● ● ● ● ● ● ● ● ● ● Z ● ● ● ● ● ●● ● ● ● ● ● ●● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●● ●● ● ● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 0 ● ● ● ● ● ● ●● ● ● ●● ● ● ● ● ●● ● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● 0 ● ● ● ● ● ● ● ● ●● ● ●●● ● ● ●● ● ● ●●● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ●● ● ● ●● ● ● ● ● ● ●● ● ● ●●● ● ● ● ● ●● ●● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ●● ● ●● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●● ● ● ● ● ● ● ●● ● ●● ● ● ● ●● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● −1 ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●● ●● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● −1 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●● ● ● ● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ●● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● −2 ● ● ● ● ● ● ● ● ● ● −2 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● −3 ● ● −3 ● −3 −2 −1 0 1 2 3 4 −3 −2 −1 0 1 2 3 4 X X Slika 1.4: Razsevni diagram realizacij X in Y ter X in Z. d) Ali je vsota X + Z porazdeljena normalno? 20 Verjetnost in statistika 1.5. Vsota normalnih spremenljivk Slika 1.5 jasno kaže, da porazdelitev vsote ni normalna. Vsota dveh normalnih spremenljivk torej ni nujno normalna, če spremenljivki nista neodvisni (dana naloga je protiprimer). 200 150 enca 100 Frekv 50 0 −4 −2 0 2 4 6 X+Z Slika 1.5: Porazdelitev vsote X + Z. Povzetek • V nalogi smo najprej izpeljali, da je vsota dveh neodvisnih standardno normalno porazdeljenih spremenljivk zopet normalna spremenljivka. Pri dokazu smo se računanju zapletenega integrala izognili z uporabo dejstva, da je integral gostote neke spremenljivke na celotni realni osi vedno enak 1 ne glede na obliko porazdelitve oziroma vrednost njenih parametrov. Morali smo torej člene preurediti tako, da je pod integra- lom ostala le gostota neke porazdelitve. • Rezultat, ki smo ga za dve spremenljivki izpeljali v nalogi, je v resnici še dosti bolj splošen - vsota poljubnega števila neodvisnih normalno porazdeljenih spremenljivk je zopet normalno porazdeljena spremen- ljivka. • Naloga je s protiprimerom pokazala, da vsota odvisnih normalno po- razdeljenih spremenljivk ni nujno normalna. Definirali smo slučajni spremenljivki X in Z in za obe pokazali, da sta normalni ter med se- boj odvisni. Njuna vsota očitno ni normalno porazdeljena. V splošnem lahko o porazdelitvi vsote odvisnih spremenljivk povemo le malo - način odvisnosti bo ključno vplival na obliko dobljene porazdelitve. 21 Verjetnost in statistika 1.6. Porazdelitev vzorčnega povprečja 1.6 Porazdelitev vzorčnega povprečja Vrnimo se spet k primeru odkrivanja dopinga. Izkaže se, da ima vsak posame- znik sebi lastno povprečje hemoglobina in da se te vrednosti med posamezniki precej razlikujejo. Da bi dosegli večjo občutljivost testa, uvedemo polletno testno obdobje, v katerem vsakega športnika testiramo petkrat. Povprečje teh petih meritev bomo vzeli kot oceno za posameznikovo povprečje pri tes- tih v prihodnosti (meje bomo postavljali glede na to povprečje). Recimo, da vemo, da se vrednosti vsakega športnika okrog njemu lastnega povprečja porazdeljujejo normalno z varianco σ2 = 52 in da so posamezne vrednosti med seboj neodvisne. a) Naj bodo Xi, i = 1, . . . , n, neodvisne, enako porazdeljene slučajne spremenljivke. Kaj lahko rečete o pričakovani vrednosti in varianci njihovega povprečja? Označite E(Xi) = µ in var(Xi) = σ2 za vsak i. n Naj bo ¯ X = 1 X n i: i=1 n n 1 1 E[ ¯ X] = E[ Xi] = E[Xi] = µ n n i=1 i=1 n n 1 1 σ2 var[ ¯ X] = var[ Xi] = var(Xi) = n n2 n i=1 i=1 Opomba: neodvisnost smo potrebovali pri izračunu variance, medtem ko bo pričakovana vrednost vsote vedno vsota pričakovanih vrednosti. b) Izračunajte meje okrog ocenjenega povprečja, znotraj katerih naj bi pri šesti meritvi nedopingiran športnik ostal z verjetnostjo 0,99 (prvih pet meritev uporabimo za oceno športnikovega povprečja). Namig: uporabite rezultat, da je vsota neodvisnih normalno porazdelje- nih spremenljivk spet normalna. Vrednosti posameznika so porazdeljene kot X ∼ N (µ, σ2) (σ = 5). Vsota 5 X i=1 i je normalna spremenljivka, prav tako je normalno porazdeljeno tudi povprečje ¯ X, saj vsoto le pomnožimo s konstanto. Zanima nas od- stopanje šeste meritve od ocenjenega povprečja v prvih petih meritvah, torej razlika Z = X 5 6 − 1 X 5 i=1 i. To je torej razlika dveh normalnih 22 Verjetnost in statistika 1.6. Porazdelitev vzorčnega povprečja spremenljivk z enakim povprečjem, ena ima varianco σ2, druga pa σ2/n. Spremenljivka Z je porazdeljena kot Z ∼ N (0, σ2 + σ2/n) = N (0,30). √ Vrednost z 5 0,005 = 2,57, meje so torej 1 X 30. 5 i=1 i ± 2,57 · c) Ali je povprečje neodvisnih enako porazdeljenih slučajnih spremenljivk vedno porazdeljeno z isto porazdelitvijo kot vsaka posamezna spremen- ljivka? Ne, to v splošnem ni res. Protiprimer je na primer vsota Bernoullijevo porazdeljenih spremenljivk, ki je porazdeljena binomsko, in prav tako to očitno ne velja za povprečje enakomerno porazdeljenih spremenljivk. Glej tudi nalogo 1.9. Predlogi za vaje v R • Predpostavite, da so povprečja športnikov normalno porazdeljena z N (148, 7,52) in generirajte povprečne vrednosti za 100 športnikov. Nato uporabite normalno porazdelitev N (0, 52), ki predstavlja odmike od osebnega povprečja vsakega športnika - generirajte po 6 vrednosti na posameznika. Ocenite osebna povprečja s pomočjo prvih petih vred- nosti ter primerjajte varianco teh povprečij s teoretično vrednostjo. Oglejte si porazdelitev odstopanja šeste vrednosti od prvih petih. Povzetek • V tej nalogi se prvič srečamo s pojmoma populacije in vzorca, s ka- terima opišemo cilj statističnega sklepanja: s pomočjo vrednosti na vzorcu želimo opisati populacijo. Pri tem si vzorec zamišljamo kot ne- odvisne realizacije neke slučajne spremenljivke, ki predstavlja popula- cijo. Vzorčenju se bomo sicer povsem posvetili v naslednjem poglavju. • Vzorčno povprečje je slučajna spremenljivka - na vsakem vzorcu iz neke populacije lahko pričakujemo drugačno vrednost. Cilj naloge je bil najti porazdelitev vzorčnega povprečja, če vemo, da je populacija normalno porazdeljena. Vzorčno povprečje je utežena vsota neodvisnih normal- nih slučajnih spremenljivk, in zato zopet normalno porazdeljena spre- menljivka. Njena pričakovana vrednost je kar populacijsko povprečje µ, njena varianca pa je enaka σ2/n. • Varianco vzorčnega povprečja imenujemo standardna napaka - varia- bilnost vzorčnega povprečja okrog populacijske vrednosti nam namreč 23 Verjetnost in statistika 1.7. Pogojna pričakovana vrednost in varianca govori o tem, koliko se bomo motili, če bomo z vzorčnim povprečjem skušali ocenjevati povprečje v populaciji. • Velikost standardne napake je obratno sorazmerna z velikostjo vzorca - večji vzorec imamo, manj se motimo pri oceni populacijskega povprečja. Standardna napaka je odvisna tudi od vrednosti σ - če je variabilnost v populaciji majhna, se pri oceni vzorčnega povprečja ne moremo dosti zmotiti. 1.7 Pogojna pričakovana vrednost in varianca Raziskovalci na področju športa so dokazali, da je pri kolesarjih hemoglobin izven tekmovalnega obdobja porazdeljen kot N (150, 72), med tekmovalnim obdobjem pa kot N (140, 112). Vzemimo, da tekmovalno obdobje traja 9 me- secev. Zanimata nas pričakovana vrednost in standardni odklon naključnega odvzetega vzorca. Namig: zanima nas slučajna spremenljivka Y , vemo {Y |X = 0} ∼ N (150, 72) in {Y |X = 1} ∼ N (140, 112), P (X = 1) = 0,75 a) Skicirajte porazdelitev Y , kaj lahko rečete o pričakovani vrednosti ter standardnem odklonu? Porazdelitev je na sliki 1.6. Pričakujemo, da bo povprečna vrednost ne- kje med vrednostima v obeh obdobjih. Ker je verjetnost, da je meritev izšla iz tekmovalnega obdobja, večja, pričakujemo, da bo vrednost ne- koliko bližja povprečju v tem obdobju. Varianca bo nekoliko večja od posameznih varianc, odvisna bo tudi od razdalje med povprečjema. b) Na danem primeru razložite formulo E(Y ) = E[E(Y |X)]. Je E(Y |X) slučajna spremenljivka ali konstanta? Izračunajte pričakovano vrednost spremenljivke Y . Z = E(Y |X) je slučajna spremenljivka, ki lahko zavzame dve vredno- sti: P (Z = 140) = 0,75, P (Z = 150) = 0,25. Pričakovana vrednost te spremenljivke je torej E(Z) = 140·P (Z = 140)+150·P (Z = 150) = 140·0,75+150·0,25 = 142,5 24 Verjetnost in statistika 1.7. Pogojna pričakovana vrednost in varianca 150000 enca 100000 Frekv 50000 0 100 120 140 160 180 Y Slika 1.6: Porazdelitev nove spremenljivke Y . Torej E[E(Y |X)] = E(Y |X = x) · P (X = x) x c) Izračunajte varianco Y . Pri izračunu variance si bomo pomagali s formulo var(Y ) = var[E(Y |X)] + E[var(Y |X)] V prvem delu gornje formule nas torej zanima varianca slučajne spremen- ljivke Z = E(Y |X): var(Z) = E([Z − E(Z)]2) = 7,52 · P [(Z − E(Z)) = 7,5] + 2,52 · P [(Z − E(Z)) = 2,5] = 56,25 · 0,25 + 6,25 · 0,75 = 18,75 Standardni odklon povprečij v različnih obdobjih okrog robnega povprečja je torej 4,33. Člen E[var(Y |X)] je pričakovana vrednost za varianco Y pri znanem X. Vemo, da je var(Y |X = 0) = 49 in var(Y |X = 1) = 121. Pričakovana vrednost je E[var(Y |X)] = 49 · 0,25 + 121 · 0,75 = 103 25 Verjetnost in statistika 1.7. Pogojna pričakovana vrednost in varianca Sestavimo oba dela skupaj in dobimo var(Y ) = 121,75, sd(Y ) = 11,03. Vrednosti izven tekmovalnega obdobja torej le malo povečajo variabilnost rezultatov. d) Izrazite varianco v splošnem: {Y |X = 0} ∼ N (µ0, σ2), {Y |X = 1} ∼ 0 N (µ1, σ2), P (X = 1) = p. 1 Vrednost E(Y |X = 0) = µ0 je pričakovana vrednost Y pri X = 0, to- rej izven tekmovalnega obdobja, podobno z µ1 = E(Y |X = 1) označimo pričakovano vrednost Y med tekmovalnim obdobjem. Verjetnost, da je športnik v tekmovalnem obdobju, označimo s p. Ker je X Bernoullijevo porazdeljena spremenljivka, velja E(X) = P (X = 1) = p. Funkcijo µ E(Y |X) = 0; X = 0 µ1; X = 1 zapišemo kot E(Y |X) = µ0(1 − X) + µ1X. Pričakovana vrednost slučajne spremenljivke Z = E(Y |X) je E(Z) = E(E(Y |X)) = E(Y |X = x)P (X = x) = µ0(1 − p) + µ1p X=x Pričakovana vrednost Y je torej uteženo povprečje pogojnih pričakovanih vrednosti v posameznih obdobjih, bližja je tistemu obdobju, iz katerega smo bolj verjetno dobili meritev. Varianca slučajne spremenljivke Z je enaka var(Z) = [E(Y |X = x) − E(Y )]2 P (X = x) X=x = [µ0 − µ0(1 − p) − µ1p]2(1 − p) + [µ1 − µ0(1 − p) − µ1p]2p = [−p(µ0 − µ1)]2(1 − p) + [(1 − p)(µ1 − µ0)]2p = [µ1 − µ0]2p2(1 − p) + [µ1 − µ0]2(1 − p)2p = [µ1 − µ0]2p(1 − p)(p + 1 − p) = [µ1 − µ0]2p(1 − p) Izrazimo še drugi del, slučajna spremenljivka var(Y |X) je enaka σ2; z verjetnostjo (1 − p) var(Y |X) = 0 σ2; z verjetnostjo p 1 26 Verjetnost in statistika 1.7. Pogojna pričakovana vrednost in varianca Spremenljivka var(Y |X) je torej Bernoullijevo porazdeljena, njena priča- kovana vrednost je E(var(Y |X)) = σ2(1 − p) + σ2p. 0 1 Združimo oba dela skupaj in dobimo var(Y ) = [µ1 − µ0]2p(1 − p) + σ2(1 − p) + σ2p 0 1 Varianca slučajne spremenljivke Y je torej seštevek uteženega povprečja varianc v posameznih obdobjih in člena, ki je odvisen od razmika med povprečjema. e) Izračunajte kovarianco X in Y . Izrazite jo splošno ({Y |X = 0} ∼ N (µ0, σ2), {Y |X = 1} ∼ N (µ ), P (X = 1) = p). 0 1, σ2 1 Kako je kovarianca odvisna od parametrov? Kaj pa korelacija? Po definiciji je kovarianca enaka cov(X, Y ) = = E[(X − E(X))(Y − E(Y ))] = (x − E(X))(y − E(Y ))fX,Y (x, y)dxdy R R Zanima nas torej pričakovana vrednost glede na skupno porazdelitev X in Y (lahko bi pisali EX,Y ). Gostote skupne porazdelitve ne poznamo, zato poizkusimo izraz preurediti tako, da bo v njem nastopala pogojna gostota. Uporabimo torej enakost fX,Y (x, y) = fY |X(y|x)fX(x) in najprej izračunajmo integral po y cov(X, Y ) = = (x − E(X))(y − E(Y ))fY |X(y|x)dy fX(x)dx R R = (x − E(X)) (y − E(Y ))fY |X(y|x)dy fX(x)dx R R V integralu E(Y )fY |X(y|x)dy lahko vrednost E(Y ) izpostavimo, saj R je konstanta. Funkcija fY |X(y|x) predstavlja pogojno gostoto - pri vsaki vrednosti x imamo torej neko slučajno spremenljivko U = Y |X=x z go- stoto fU (u) = fY |X(y|x). Zato je integral pri dani vrednosti x enak 27 Verjetnost in statistika 1.7. Pogojna pričakovana vrednost in varianca fY |X(y|x)dy = 1, torej R cov(X, Y ) = = (x − E(X)) yfY |X(y)dy − E(Y ) fX(x)dx R R = (x − E(X)) [E(Y |X) − E(Y )] fX(x)dx R V našem primeru je X diskretna spremenljivka, integral po X lahko torej zamenjamo z vsoto dveh členov: cov(X, Y ) = = (0 − E(X))[E(Y |X = 0) − E(Y )]P (X = 0) + (1 − E(X))[E(Y |X = 1) − E(Y )]P (X = 1) Vrednost E(Y |X = 0) = µ0 je pričakovana vrednost Y pri X = 0, torej izven tekmovalnega obdobja, podobno smo z µ1 = E(Y |X = 1) označili pričakovano vrednost Y med tekmovalnim obdobjem. Označimo še E(Y ) = µ. Verjetnost, da je športnik v tekmovalnem obdobju, označimo s p. Ker je X Bernoullijevo porazdeljena spremenljivka, velja E(X) = P (X = 1) = p. Zato cov(X, Y ) = = −p[µ0 − µ](1 − p) + (1 − p)[µ1 − µ]p = p(1 − p)(−µ0 + µ1) in p(1 − p)(µ cor(X, Y ) = 1 − µ0) √varY p(1 − p) V našem primeru cov(X, Y ) = 0,75 · 0,25 · (140 − 150) = −1,875 1,875 cor(X, Y ) = − √ = −0,392 11,03 · 0,75 · 0,25 Kovarianca in korelacija sta odvisni od razlike med povprečjema - večja kot je razlika, večji sta (po absolutni vrednosti). Če bi bila razlika 0, torej 28 Verjetnost in statistika 1.7. Pogojna pričakovana vrednost in varianca povprečje neodvisno od obdobja, bi bili tudi korelacija oz. kovarianca enaki 0. Vrednosti sta negativni, kadar večji X pomeni manjši Y . Odvisni sta tudi od p - največji sta, kadar je p = 0,5, torej kadar obe obdobji enako prispevata k skupnemu povprečju (če bi bilo vrednosti v enem obdobju zelo malo, bi bila korelacija majhna). Dodatno je korelacija odvisna tudi od variabilnosti v enem in drugem obdobju. Če bi bila ta variabilnost velika v primerjavi z razliko med povprečjema, spremenljivki ne bi bili močno povezani. f) Kolikšne so vrednosti variance, kovariance in korelacije, če sta povprečji v tekmovalnem in izven tekmovalnega obdobja enaki? Ali sta spremenljivki X in Y tedaj neodvisni? Če je razlika enaka 0, torej povprečje neodvisno od obdobja, je varianca Y enaka var(Y ) = σ2(1 − p) + σ2p, korelacija in kovarianca pa sta enaki 0 1 0. Vendar pa to ne pomeni, da sta spremenljivki X in Y neodvisni - od vrednosti X je odvisna varianca Y . Porazdelitev Y je torej odvisna od X, četudi X ne vpliva na povprečje. Torej: vemo, da je korelacija enaka 0, če sta spremenljivki neodvisni, vendar obratno ni nujno res. Predlogi za vaje v R • Generirajte veliko število vrednosti in si oglejte njihovo porazdelitev: > set.seed(1) > a <- rnorm(1000000,mean=140,sd=11) #vrednosti Y|X za tek. obd. > b <- rnorm(1000000,mean=150,sd=7) #vrednosti Y|X za ne-tek. obd. > x <- sample(0:1,size=1000000, #obdobje - porazdelitev X + replace=T,prob=c(0.25,0.75)) > y <- a*x+b*(1-x) #slucajna spremenljivka Y > hist(y,prob=T) #narisemo spremenljivko > mean(y) #ocena povprecja > var(y) #ocena variance • Poizkusite preveriti vsakega od rezultatov še z R. Primerjajte teoretične vrednosti z njihovimi ocenami. Povzetek • Pogojna pričakovana vrednost E(Y |X) je slučajna spremenljivka. Pri- čakovana vrednost te slučajne spremenljivke je E[E(Y |X)] = E(Y ). 29 Verjetnost in statistika 1.8. Uvrščanje • Če poznamo pogojno porazdelitev Y |X, lahko varianco spremenljivke Y izrazimo s formulo var(Y ) = var[E(Y |X)] + E[var(Y |X)] • Kovarianca in korelacija dveh neodvisnih spremenljivk sta enaki 0, obratno pa ni nujno res: četudi je kovarianca oziroma korelacija enaka 0, spremenljivki nista nujno neodvisni. Protiprimer je razložen v točki f). 1.8 Uvrščanje Na podlagi nekega kazalnika želimo ocenjevati kreditno sposobnost posamez- nikov, želimo jih uvrstiti v dve skupini - tiste, ki bodo kredit odplačali, in tiste, ki ga ne bodo. Kot učni vzorec imamo na razpolago vrednosti tega kazalnika za posameznike, ki so lansko leto najeli enoletni kredit, in podatke o tem, ali je letos kredit odplačan ali ne (naloga je povzeta po Blagus, 2011). Predpostavimo, da so vrednosti kazalnika porazdeljene pri obeh skupinah posameznikov približno normalno. Na podlagi letošnjih podatkov ocenimo povprečno vrednost kazalnika za posameznike, ki so kredit odplačali ( ¯ Xd), in za posameznike, ki ga niso ( ¯ Xs). Ti dve oceni sedaj uporabimo za uvrščanje strank: napovemo, da bo nek posameznik uspel odplačati kredit, če je tre- nutna vrednost njegovega kazalnika bližja ¯ Xd kot ¯ Xs. Raziskati želimo lastnosti takega uvrščanja. Recimo, da smo v učni vzorec zajeli nd = 750 posameznikov, ki so kredit odplačali, in ns = 250, ki ga niso. Recimo, da je kazalnik povsem neuporaben za naš namen, torej da je nje- gova porazdelitev enaka pri “dobrih” kot pri “slabih” strankah: Xs in Xd sta enako porazdeljena, označimo kar z X: X ∼ N (50, 152). Ugotoviti želimo, kolikšna bo verjetnost, da neko naključno stranko na podlagi današnje vred- nosti njenega kazalnika uvrstimo med “dobre”. a) Kakšna je porazdelitev ¯ Xs (v splošnem, torej za neko varianco σ2 in neko število slabih ns v učnem vzorcu)? Xs ∼ N (µ, σ2). Vemo, da je povprečje neodvisnih normalnih spremenljivk normalno porazdeljeno ter da velja ¯ Xs = 1 X ). n s,i ∼ N (µ, σ2 s ns 30 Verjetnost in statistika 1.8. Uvrščanje b) Kakšna je porazdelitev X − ¯ Xs? Vemo že (glej nalogo 1.6), da velja X − ¯ Xs ∼ N (0, σ2 + σ2 ). ns c) Označimo Z = X − ¯ Xs in Y = X − ¯ Xd. Izračunajte kovarianco in kore- lacijo. Interpretirajte izračunano formulo za korelacijo za ns = nd. Kako je odvisna od velikosti vzorcev? Nova vrednost X je neodvisna od vzorčnih povprečij, povprečji pa sta med seboj prav tako neodvisni. Zato velja cov[X − ¯ Xs, X − ¯ Xd] = cov[X, X] = var(X) = σ2 Korelacija je enaka σ2 cor[X − ¯ Xs, X − ¯ Xd] = σ2(1 + 1 ) σ2(1 + 1 ) ns nd 1 = (1 + 1 ) (1 + 1 ) ns nd Če sta velikosti vzorca med seboj enaki, velja n cor[X − ¯ X s s, X − ¯ Xd] = 1 + ns Ko vzorca naraščata, gre korelacija proti 1. To je intuitivno jasno, saj z naraščanjem vzorca oceni povprečij postaneta zelo natančni (v primerjavi s posamezno vrednostjo sta skoraj konstanti). d) Izrazite formulo za gostoto fZ,Y (z, y). Namig: uporabimo rezultat, da je skupna gostota v tem primeru nor- malna (Rice (2009), stran 102), kar pa v splošnem ni nujno vedno res. Za zapis dvorazsežne normalne porazdelitve potrebujemo robna povprečja in variance (izpeljali smo jih v drugi točki) ter korelacijo ρ med spremen- ljivkama (izpeljali smo jo v prejšnji točki). Skupno gostoto torej zapišemo 31 Verjetnost in statistika 1.8. Uvrščanje 3 ● ● ● 2 ● ● ● ● ● ● ● ● ● ● ● ● 1 ● ● ● ● ●●● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● Y 0 ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● −1 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● −2 ● ● ● −3 ● −3 −2 −1 0 1 2 3 Z Slika 1.7: Razsevni diagram razlik za 100 realizacij slučajne spremenljivke. kot 1 1 (z − µ (y − µ f Z )2 Y )2 Z,Y (z, y) = exp − + 2πσ 2(1 − ρ2) σ2 σ2 Z σY 1 − ρ2 Z Y 2ρ(z − µ − z )(y − µy ) σZσY 1 1 z2 y2 2ρzy = exp − + − 2πσ2sd 1 − ρ2 2σ2(1 − ρ2) s2 d2 sd Pri tem smo označili s = 1 + 1 in d = 1 + 1 , torej ρ2 = 1/(sd). ns nd Definirajmo še c = nsnd in dobimo 1+ns+nd c − c2 (y2 1+ns +z2 1+nd −2yz) f 2σ2 ns n Z,Y (z, y) = e d 2πσ2 e) Zanima nas verjetnost P (|X − ¯ Xs| < |X − ¯ Xd|). Kako bi jo izračunali? Skicirajte območje, ki nas zanima, nastavite integral in meje. Slika 1.8 prikazuje območje integracije. Integral, ki ga moramo izračunati, je enak P (|Y | > |Z|) = fZ,Y (z, x)dxdz |Y |>|Z| c − c2 (y2 1+ns +z2 1+nd −2yz) = e 2σ2 ns nd dzdy 2πσ2 |Y |>|Z| 32 Verjetnost in statistika 1.8. Uvrščanje Y3 ● ● ● ● 2 ● ● ● ● ● ● ● ● ● ● ● ● 1 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● 0 ●● ● ● ● ● ● ●● ● −3 −2 ● −1 0 ●● 1 2 3 Z ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● −1 ● ● ● ● ● ● ● −2 ● ● ● −3 Slika 1.8: Območje |Z| < |Y |, za ns = 10, nd = 20. Zaradi simetrije je integral za pozitivne y (glej sliko 1.8) enak kot za negativne, zato je dovolj, da integriramo le gornji del (in množimo z 2). Gornja kraka dodatno razdelimo na negativne in pozitivne vrednosti z: c ∞ y − c2 (y2 1+ns +z2 1+nd −2yz) P (|Y | > |Z|) = e 2σ2 ns nd dzdy + πσ2 0 0 c ∞ 0 − c2 (y2 1+ns +z2 1+nd −2yz) + e 2σ2 ns nd dzdy. πσ2 0 −y Izračunamo zgornje integrale in dobimo 1 1 n P (|Y | > |Z|) = arctan dns + π nd 1 + nd + ns 1 + 2n n + arctan d dns . nd 1 + nd + ns Za ns = nd je rezultat 0,5, za ns = 250, nd = 750 pa 0,49. Če sta vzorca enako velika, bo ta način uvrščanja enote v vsako skupino uvrstil z ver- jetnostjo 0,5, kar je smiselno (saj kazalnik nič ne pove o kvaliteti). Čim pa vzorca nista enako velika, verjetnost ne bo več enaka 0,5 niti ne bo proporcionalna velikosti vzorca. Za neenake velikosti vzorcev (neurav- notežene podatke), to uvrščanje torej ne daje želenih oz. pričakovanih rezultatov (je pa res, da bodo morali biti podatki zelo neuravnoteženi, da bo verjetnost bistveno različna od 0,5). 33 Verjetnost in statistika 1.8. Uvrščanje Predlogi za vaje v R • Generirajte vrednosti obeh slučajnih spremenljivk in narišite dobljene razlike. Pomagajte si s spodnjo kodo. > set.seed(1) > slabi <- rnorm(10) #10 vrednosti kazalnika za slabe > dobri <- rnorm(10) #10 vrednosti kazalnika za dobre > nov <- rnorm(1) #ena vrednost za novega posameznika > plot(nov-mean(slabi),nov-mean(dobri), + ylim=c(-3,3),xlim=c(-3,3)) #razdalja vrednosti od obeh povprecij > abline(0,1) #simetrala > for(it in 1:99){ #ponovim 100x, vsakic znova simuliram + slabi <- rnorm(10) #10 slabih in 10 dobrih + dobri <- rnorm(10) # ter enega novega, ki ga na podlagi + nov <- rnorm(1) # dobljenega kriterija uvrstim + points(nov-mean(slabi),nov-mean(dobri)) + } • Spreminjajte velikosti vzorcev in opazujte, kako se spreminja verjetnost uvrstitve v posamezno skupino. Preštejte vrednosti, ki so uvrščene v vsako skupino, in ocenite verjetnosti. • Podatke generirajte še na primeru bolj smiselnega indeksa (ki ima različni povprečji v skupinah). Preverite, kako dobro sedaj uvrščate vrednosti. Povzetek • Naloga uporabi predhodno pridobljeno znanje za izpeljavo skupne po- razdelitve dveh slučajnih spremenljivk. Verjetnost, ki nas zanima, nato izračunamo kot integral po ustreznem delu tega prostora, torej območju, kjer sta slučajni spremenljivki v želeni zvezi. • V nalogi pokažemo, da pri neuravnoteženih podatkih (vzorca različne velikosti) sicer na videz smiselno uvrščanje ne daje pričakovanih rezul- tatov. Smiselnost rezultatov smo preverili za primer, ko kazalnik ne nosi nobene informacije o skupini - ker uvrščanje že v tem primeru ni smiselno, bo rezultate uvrščanja nemogoče interpretirati v primeru, ko bo kazalnik pomemben. 34 Verjetnost in statistika 1.9. Centralni limitni izrek 1.9 Centralni limitni izrek V raziskavo o pojavih alergije pri dojenčkih so vključili vse otroke iz lokalnega zdravstvenega doma, ki so bili na nek datum stari manj kot eno leto, in jih nato spremljali pol leta. Kakšno povprečno starost otrok ob vključitvi v raz- iskavo pričakujete? Določite meje, v katerih bo s 95% verjetnostjo povprečna starost otrok na vzorcu, če so bili otroci v raziskavo zares izbrani naključno (predpostavimo, da so vsi rojstni dnevi v letu enako verjetni). a) Kakšna je približna porazdelitev vzorčnega povprečja ¯ X za nek n? Zaradi preprostosti vzemimo, da je starost X kar zvezna spremenljivka, torej enakomerno porazdeljena na intervalu [0,1]. Vsota enakomerno po- razdeljenih spremenljivk očitno ni enakomerno porazdeljena, za izpelje- vanje porazdelitve vsote bi lahko uporabili formulo za vsoto slučajnih spremenljivk (glej nalogo 1.4). Ker nam zadostuje približek, bomo ubrali lažjo pot in uporabili centralni limitni izrek, ki pravi, da bo porazdelitev vsote za dovolj velik n približno normalna. Izračunati moramo torej pa- rametra te normalne porazdelitve, ker vemo, da velja E( ¯ X) = E(X) in var( ¯ X) = varX/n, potrebujemo le pričakovano vrednost in varianco spre- menljivke X. Uporabimo, da je gostota fX(x) enakomerne porazdelitve enaka 1 na intervalu [0,1] in 0 sicer, in dobimo 1 x2 1 E(X) = xfX(x)dx = xdx = |1 = 2 0 2 R 0 1 var(X) = (x − E(X))2fX(x)dx = (x − 0,5)2dx R 0 1 1 x3 x2 x 1 = (x2 − x + )dx = ( − + )|1 = 4 3 2 4 0 12 0 Porazdelitev ¯ X torej lahko aproksimiramo z N ( 1 , 1 ). 2 12n b) Zapišite meje, v katerih bo s 95% verjetnostjo vzorčno povprečje, če v vzorec zajamemo 100 otrok. 35 Verjetnost in statistika 1.9. Centralni limitni izrek Izračunali smo, da je približna porazdelitev povprečnega časa enaka ¯ X ∼ N (0,5, 0,0008). Pri tej porazdelitvi se vrednosti spremenljivke s 95% ver- √ jetnostjo nahajajo v mejah 0,5 ± 1,96 · 0,0008, torej [0,44, 0,56]. Predlogi za vaje v R • Generirajte velik vzorec iz enakomerne porazdelitve (runif) in na njem preverite izračun variance. • Oglejte si porazdelitev vzorčnega povprečja za vzorce velikosti n = 5, 20, 50 in grafično s histogramom preverite, kako dobro se jim prilega ustrezna normalna porazdelitev. • Simulirajte vrednosti iz različnih porazdelitev (npr. rgamma, rchisq, rf, rt ...) in izračunajte njihovo povprečje. Postopek ponovite veli- kokrat in si oglejte porazdelitev vzorčnega povprečja in na ta način s simulacijami preverite centralni limitni izrek. • Poizkusite si izmisliti porazdelitev, pri kateri bo pri vzorcih velikosti n = 50 normalna porazdelitev še vedno slaba aproksimacija za poraz- delitev vzorčnega povprečja. Povzetek • Centralni limitni izrek nam poda asimptotsko aproksimacijo porazdeli- tve vsote neodvisnih enako porazdeljenih slučajnih spremenljivk. Izrek je ključnega pomena v statistični teoriji, saj pomaga izpeljati poraz- delitev neke vzorčne statistike, ne da bi morali poznati porazdelitev populacije, edini pogoj je, da se dà statistika zapisati kot vsota ne- odvisnih slučajnih realizacij iz te populacije. V tem primeru bomo porazdelitev vzorčne statistike lahko aproksimirali z normalno in to je tudi razlog, zakaj se normalna porazdelitev v statistiki tako pogosto pojavlja. • Pri uporabi centralnega limitnega izreka se moramo seveda zavedati, da gre vedno le za približek, ki bo tem boljši, čim večji bo vzorec. Kako dobra je aproksimacija na majhnem vzorcu, bi lahko preverili s simulacijami. 36 Verjetnost in statistika 1.10. Normalna aproksimacija • V nalogi 1.4 smo pokazali primer porazdelitve, pri kateri je vsota neod- visnih enako porazdeljenih slučajnih spremenljivk zopet porazdeljena z isto porazdelitvijo (spremenijo se le parametri). Po drugi strani nam centralni limitni izrek pove, da vsota neodvisnih enako porazdeljenih slučajnih spremenljivk konvergira proti normalni porazdelitvi. Z vsoto torej lahko ostajamo znotraj iste porazdelitve le pri tistih porazdeli- tvah, ki so pri določenih vrednostih parametrov zelo podobne normalni porazdelitvi oziroma konvergirajo proti normalni porazdelitvi. 1.10 Normalna aproksimacija binomske po- razdelitve V nekem kraju želijo zgraditi obvoznico, zanima jih delež ljudi, ki to gradnjo podpirajo. V ta namen izvedejo anketo. Recimo, da je verjetnost, da se posameznik strinja, enaka 0,65. Kolikšna je verjetnost, da bo med 6 posamezniki večina za gradnjo? a) Naj bo X = I{posameznik se strinja}, X je Bernoullijevo porazdeljena spremenljivka. Kako je porazdeljena vsota S6 = 6 X i=1 i? Izračunajte pričakovano vrednost in standardni odklon. Pričakovana vrednost Bernoullijeve spremenljivke je E(X) = 0 · P (X = 0) + 1 · P (X = 1) = p var(X) = E(X2) − E(X)2 = 12 · P (X = 1) − p2 = p − p2 = p(1 − p) za vsoto pa velja n n E(Sn) = E( Xi) = E(Xi) = np i=1 i=1 n n var(Sn) = var( Xi) = var(Xi) = np(1 − p) i=1 i=1 Vsota neodvisnih enako Bernoullijevo porazdeljenih spremenljivk je po- razdeljena binomsko (glej nalogo 1.3), verjetnost posameznega izida je n P (Sn = k) = pk(1 − p)n−k k 37 Verjetnost in statistika 1.10. Normalna aproksimacija b) Izračunajte verjetnost, da so vsaj 4 posamezniki za gradnjo. Z uporabo formule za binomsko porazdelitev: 6 6 P (Sn ≥ 4) = (0,65)k(0,35)6−k = 0,647 k k=4 c) Aproksimirajte to verjetnost še s pomočjo centralnega limitnega izreka. Vemo, da Sn−np √ konvergira (v porazdelitvi) proti standardno normalno np(1−p) porazdeljeni spremenljvki Z. Poglejmo, kako dobra bo aproksimacija z normalno porazdelitvijo pri n = 6: S 3,5 − np P (S n − np n > 3,5) = P ( > ) np(1 − p) np(1 − p) 3,5 − 3,9 = P (Z > ) = 0,634 1,17 Na slikah 1.9 in 1.10 je prikazana kvaliteta aproksimacije za dve vrednosti p. Vidimo, da je aproksimacija nekoliko slabša, če je porazdelitev bolj asimetrična, a še vedno zelo dobra. 0.4 ● 0.10 0.3 0.08 ● ● 0.2 0.06 0.04 0.1 ● ● 0.02 ● ● 0.0 0.00 Slika 1.9: Aproksimacija binomske porazdelitve z normalno za p = 0,65 in (a)n = 6, (b)n = 500. 38 Verjetnost in statistika 1.10. Normalna aproksimacija 0.6 0.5 ● 0.15 0.4 ● 0.10 0.3 0.2 0.05 0.1 ● 0.0 ● ● ● ● 0.00 Slika 1.10: Aproksimacija binomske porazdelitve z normalno za p = 0,9 in (a)n = 6, (b)n = 500. Predlogi za vaje v R • Oglejte si, kako dobra je aproksimacija za različne velikosti vzorca in različne vrednosti p: > win.graph(height=6,width=12) #pripravimo okno za risanje > par(mfrow=c(1,2)) #narisali bomo dva grafa na isto sliko > p <- 0.65 #verjetnost, da je posameznik za > n <- 6 #velikost vzorca > verb <- dbinom(0:n,n,p) #verjetnosti posameznih izidov > barplot(verb,space=0,width=1,ylim=c(0,.4)) #stolpicni diagram > sd <- sqrt(n*p*(1-p)) #standardni odklon > e <- n*p #pricakovana vrednost > vern <- dnorm(0:n,e,sd) #vrednost gostote v posameznih tockah > points(0:n+.5,vern,pch=16) #dorisemo vrednosti na graf > sek <- seq(0,n+1,length=100) #izberemo tocke za racunanje gostote > vern <- dnorm(sek,e,sd) #vrednost gostote v izbranih tockah > lines(sek+.5,vern,lwd=2) #dorisemo krivuljo na graf ################## > n <- 50 #se enkrat za vecji vzorec > verb <- dbinom(0:n,n,p) #verjetnosti posameznih izidov > barplot(verb,space=0,width=1) #stolpicni diagram > sd <- sqrt(n*p*(1-p)) #standardni odklon > e <- n*p #pricakovana vrednost > vern <- dnorm(0:n,e,sd) #vrednost gostote v posameznih tockah > lines(0:n+.5,vern,lwd=2) #dorisemo vrednosti na graf Povzetek 39 Verjetnost in statistika 1.10. Normalna aproksimacija • Centralni limitni izrek nam poda asimptotsko aproksimacijo porazdeli- tve neke spremenljivke, ki se jo dà zapisati kot vsoto neodvisnih enako porazdeljenih slučajnih spremenljivk. V tej nalogi smo porazdelitev vsote poznali in smo izrek uporabili le kot pripomoček (približek) za lažje oziroma hitrejše računanje. • Da računanje z binomsko porazdelitvijo ne bi bilo prezamudno, smo uporabili zelo majhen n. Navkljub temu je bil približek z normalno porazdelitvijo že kar soliden, seveda bi bil za večje vrednosti n še toliko boljši. 40 Poglavje 2 Vzorčenje Statistični del knjige pričenjamo s poglavjem o vzorčenju (Rice, 2009, po- glavje 7), ki bo služilo kot osnova za razumevanje lastnosti cenilk. Z izrazom cenilka označujemo funkcijo slučajnih spremenljivk (vrednosti na vzorcu), s katero želimo oceniti populacijsko količino, ki nas zanima. Pri vsaki cenilki nas bo zanimalo predvsem dvoje - ali je cenilka nepristranska in kolikšna je njena varianca. Poiskali bomo nepristranske cenilke za različne primere po- pulacij in načine vzorčenja - začeli bomo z neskončno populacijo, saj slučajno vzorčenje iz take populacije da neodvisne enako porazdeljene slučajne spremenljivke. Nato si bomo ogledali, kako se izrazi za cenilke in njihove variance spremenijo, če se med spremenljivkami pojavi odvisnost bodisi zaradi končne populacije bodisi zaradi dejstva, da so enote vzorčene iz posameznih skupin s specifičnimi lastnostmi. V večini nalog se bomo osredotočili na ocenjevanje povprečja - proučevali bomo lastnosti intuitivno smiselnih cenilk populacijskega povprečja. Naša prva naloga bo pokazati, da je cenilka nepristranska, torej da je njena pričakovana vrednost enaka povprečju v populaciji. Naslednji korak je poiskati standardno napako te ocene, skušali jo bomo izraziti z varianco populacije ter velikostjo vzorca. Tako izraženo varianco lahko potem primerjamo med alternativnimi cenilkami in ugotavljamo, katera ima najmanjšo. Zadnji ko- rak večine nalog je poiskati še cenilko variance, torej ugotoviti, kako njeno velikost lahko ocenimo s pomočjo vzorca. 41 Verjetnost in statistika 2.1. Vzorčenje - neskončna populacija 2.1 Vzorčenje - neskončna populacija Oceniti želimo znižanje vrednosti pritiska pri pacientih z esencialno hipertenzijo po treh mesecih jemanja nekega zdravila. V ta namen smo zbrali vzorec 25 bolnikov, naj bo Xi vrednost razlike pri i-tem bolniku našega vzorca. Predpostavimo, da so slučajne spremenljivke Xi neodvisne in enako poraz- deljene. a) Pokažite, da je povprečje našega naključnega vzorca nepristranska ocena povprečnega znižanja v populaciji bolnikov (to označimo z µ). Povprečje vzorca je 1 n X n i=1 i, povprečje populacije označimo z µ. Pred- postavljamo, da je vzorec naključen, torej da so vrednosti Xi enako po- razdeljene, za vse velja E(Xi) = µ. n n 1 1 1 E( Xi) = E(Xi) = nµ = µ n n n i=1 i=1 b) Kaj lahko rečemo o cov(Xi, Xj), če je i = j? Ker so vrednosti bolnikov med seboj neodvisne, je kovarianca enaka 0. c) Naj bo varianca v populaciji enaka σ2. Kolikšna je varianca (oz. standar- dna napaka) naše cenilke? n n n n 1 1 1 var ¯ X = cov Xi, Xj = cov Xi, Xj n n n2 i=1 j=1 i=1 j=1 n n 1 = cov Xi, Xj n2 i=1 j=1 n n 1 = cov (Xi, Xi) + cov (Xi, Xj) n2 i=1 j=1,j=i n 1 = [cov (Xi, Xi) + (n − 1) cov (Xi, Xj)] n2 i=1 42 Verjetnost in statistika 2.1. Vzorčenje - neskončna populacija Uporabimo, da so vrednosti med seboj neodvisne, torej da je cov[Xi, Xj] = 0 za vsak i = j. n 1 var ¯ X = [cov(Xi, Xi)] n2 i=1 1 1 = · ncov(Xi, Xi) = var(Xi) n2 n σ2 = n Standardna napaka je enaka SE = σ √ . n d) Na podlagi našega vzorca želimo oceniti σ2. Naj bo naša cenilka σ2 = c n (X i=1 i − ¯ X)2. Kolikšna mora biti vrednost konstante c, da bo naša cenilka nepristranska? Vemo, da velja σ2 = E(X2) − E(X)2, torej E(X2) = σ2 + µ2 in podobno tudi SE2 = var( ¯ X) = σ2 = E( ¯ X2) − E( ¯ X)2, torej E( ¯ X2) = µ2 + σ2 . n n Pri prehodu iz prve v drugo vrstico upoštevamo, da velja X ¯ X = i i ¯ X X ¯ X2: i i = n ¯ X2 = i n n 2 E c X ¯ i − ¯ X = cE X2 − 2X X + ¯ X2 i i i=1 i=1 n = cE X2 − n ¯ X2 i i=1 n = cE X2 − ¯ X2 i i=1 n = c E(X2) − E( ¯ X2) i i=1 n σ2 = c µ2 + σ2 − µ2 + n i=1 1 = cn σ2 1 − n = σ2(n − 1)c 43 Verjetnost in statistika 2.1. Vzorčenje - neskončna populacija Ker želimo, da velja E(σ2) = σ2, mora biti c = 1 . n−1 e) Na vzorcu smo dobili naslednje rezultate: ¯ x = 4, ˆ σ = 20. Ocenite stan- dardno napako (torej standardni odklon vzorčnega povprečja). Ali boste na podlagi podatkov lahko trdili, da se tlak zniža tudi v populaciji? Ocena standardne napake je enaka SE = ˆσ √ = 20 = 4. Znižanje pri- n 5 tiska je enako standardni napaki - prave razlike ne moremo razločiti od naključne variabilnosti. V ta namen bi potrebovali precej večji vzorec. f) Kako bi s pomočjo dobljenih podatkov zapisali 95% interval zaupanja za populacijsko povprečje? Interval zaupanja je verjetnostni interval, za njegov izračun moramo po- znati porazdelitev cenilke. Če bi bile vrednosti pritiska v populaciji nor- malno porazdeljene, bi bilo tudi vzorčno povprečje normalno porazdeljeno (naloga 1.6). Ker smo varianco cenilke ocenili iz podatkov, bi za meje in- tervala zaupanja morali uporabiti porazdelitev t (Rice, 2009, poglavje 6). Če pa porazdelitev v populaciji ni znana, si pomagamo s centralnim li- mitnim izrekom, ki nam pove, da lahko porazdelitev vzorčnega povprečja aproksimiramo z normalno. Dobimo torej (glej npr. Rice, 2009, razdelek 7.3.3): ¯ X ± zα/2 · SE = 4 ± 1,96 · 4 = [−3,84, 11,84] Predlogi za vaje v R • Generirajte vzorce velikosti 30 iz enakomerne porazdelitve, vsakič iz- računajte povprečje ter si oglejte porazdelitev tega povprečja. Ocenite pričakovano vrednost za povprečje in standardno napako ter oceni pri- merjajte s teoretičnimi vrednostmi. • Postopek ponovite še za druge (morda bolj asimetrične) porazdelitve. Povzetek • Medsebojna neodvisnost enot v primeru neskončne populacije močno pripomore k preprostosti izpeljav. • Nepristransko cenilko za varianco smo izpeljali tako, da smo najprej zapisali intuitivno smiseln izraz, nato pa določili še konstanto, pri kateri je ta cenilka nepristranska. 44 Verjetnost in statistika 2.2. Vzorčenje - končna populacija • Če je populacija neskončna in enote v vzorcu neodvisne, velja: n σ2 1 var( ¯ X) = ; σ2 = (Xi − ¯ X)2 n n − 1 i=1 • Da bi na podlagi vzorca lahko napovedovali vrednost povprečja v popu- laciji, bi seveda morali poznati celotno porazdelitev in ne le pričakovane vrednosti ter variance vzorčnega povprečja. Če so vrednosti X nor- malno porazdeljene, je vzorčno povprečje prav tako normalno porazde- ljeno, da bi lahko računali poljubne verjetnosti, moramo zato izračunati le parametra te porazdelitve. A tudi v splošnem lahko porazdeli- tev vzorčnega povprečja aproksimiramo z normalno porazdelitvijo - tu ključno vlogo odigra centralni limitni izrek. 2.2 Vzorčenje - končna populacija Oceniti želimo povprečno število zaposlenih v podjetjih neke panoge ob začet- ku letošnjega leta. Panoga je razdeljena na podskupine, v eni izmed skupin je le 11 podjetij. Uspeli smo pridobiti podatke za naključen vzorec šestih izmed teh podjetij. Naj bo Xi število zaposlenih v i-tem podjetju našega vzorca, µ naj označuje njihovo povprečje, σ pa standardni odklon. a) Naj X1 in X2 označujeta vrednosti prvih dveh naključno izbranih podje- tij. Kaj lahko rečemo o cov(X1, X2)? Kaj pa za splošen i = j? Populacija je končna, označimo vrednosti v populaciji z xk, k = 1, . . . ,11. Po nekem vrstnem redu izberimo vseh 11 podjetij, prvih 6 naj predsta- vlja vzorec, Xi označuje število zaposlenih v i-tem izbranem podjetju. Ker je vsak izmed Xi ena od vrednosti xk in imajo vse enako verjetnost, je cov(X1, X2) = cov(Xi, Xj) za poljubna različna i in j. Vendar pa sedaj Xi in Xj nista neodvisni slučajni spremenljivki - če je Xi = xk, Xj ne more zavzeti k-te vrednosti. N N cov Xi, Xj = cov Xi, xk = 0 j=1 k=1 45 Verjetnost in statistika 2.2. Vzorčenje - končna populacija Ker je vsota vseh vrednosti konstanta, je zgornji izraz enak 0, torej N cov Xi, Xj = cov (Xi, Xi) + (N − 1)cov (Xi, Xj) = 0 j=1 in zato (za i = j) σ2 cov(Xi, Xj) = − N − 1 Spremenljivki sta torej negativno povezani. b) Izračunajte še korelacijo cor(Xi, Xj) za neka i = j. Od česa je odvisna? Korelacija med spremenljivkama je enaka σ2 1 cor(Xi, Xj) = − = − (N − 1)σ2 N − 1 Korelacija je torej odvisna izključno od velikosti populacije in z večanjem populacije hitro postane zelo majhna. c) Izračunajte standardno napako vzorca velikosti n. n n n n 1 1 1 1 var[ ¯ X] = cov[ Xi, Xi] = ·cov[Xi, Xi] n n n2 n i=1 i=1 i=1 i=1 n 1 = {cov[Xi, Xi] + (n − 1)cov[Xi, Xj]} n2 i=1 n 1 σ2 σ2 N − n = {σ2 − (n − 1) } = n2 N − 1 n (N − 1) i=1 Standardna napaka je torej enaka SE = σ √ N −n . n N −1 d) V drugi panogi imamo 100 podjetij. Kako velik vzorec moramo vzeti iz te panoge, da bomo dobili približno enako veliko standardno napako (pri- vzemimo, da je varianca tudi v tej panogi enaka σ2)? Kaj pa če bi imeli panogo z zelo velikim številom podjetij? 46 Verjetnost in statistika 2.2. Vzorčenje - končna populacija Za N = 11 in n = 6 dobimo SE2 = σ2 11−6 = σ2 . Pri populaciji ve- 6 (10) 12 likosti 100 nam vzorec velikosti 10 da standardno napako SE2 = σ2 , 11 vzorec velikosti 11 pa standardno napako SE2 = σ2 . 12,2 Če je populacija večja, bomo za enako standardno napako torej potre- bovali več enot. Če bi imeli skupino z zelo velikim številom podjetij, bi bil izraz N−n približno enak 1, zato bi potrebovali 12 podjetij. Velikost (N −1) potrebnega vzorca se z večanjem populacije veča, vendar pa to večanje kmalu ni več bistveno. e) Kolikšna mora biti vrednost konstante c, da bo cenilka σ2 = c n (X i=1 i − ¯ X)2 nepristransko ocenila vrednost σ2? Ponovimo izračun iz zadnje točke prejšnje naloge, upoštevamo, da je E( ¯ X2) = µ2 + σ2 N−n : n N −1 n n E[c (Xi − ¯ X)2] = cE X2 − ¯ X2 i i=1 i=1 n σ2 N − n = c µ2 + σ2 − µ2 + n N − 1 i=1 N − n = cn σ2(1 − ) n(N − 1) N (n − 1) = σ2c N − 1 Ker želimo, da velja E(σ2) = σ2, mora biti c = 1 N−1 , torej n−1 N n 1 N − 1 σ2 = (Xi − ¯ X)2 n − 1 N i=1 f) Zapišite še nepristransko cenilko za varianco povprečja. Združimo dosedanje rezultate in dobimo n 2 σ2(N − n) 1 N − 1 (N − n) SE = = (Xi − ¯ X)2 N − 1 n − 1 N N − 1 i=1 σ2 N − n = n N 47 Verjetnost in statistika 2.3. Določitev načrta vzorčenja kjer je σ2 = 1 n (X n−1 i=1 i − ¯ X)2. Predlogi za vaje v R • Izmislite si 11 vrednosti ter nato generirajte vzorce velikosti 6. Oglejte si porazdelitev vzorčnih povprečij. Pokažite, da gornja formula pred- stavlja nepristransko oceno populacijske variance. Povzetek • Če je populacija končna, vrednosti enot izbranih v vzorec navkljub na- ključnemu izbiranju med seboj niso neodvisne, kovarianca med enotami je enaka: σ2 cov(Xi, Xj) = − N − 1 • Varianca vzorčnega povprečja in nepristranska cenilka za varianco v populaciji sta enaki n σ2 N − n 1 N − 1 var( ¯ X) = ; σ2 = (Xi − ¯ X)2 n (N − 1) n − 1 N i=1 • V končni populaciji bo na izbiro velikosti vzorca vplivala tudi veli- kost populacije. Vendar pa bo ta vpliv pomemben le pri sorazmerno majhnih populacijah oziroma takrat, ko velikost vzorca predstavlja ne- zanemarljiv delež velikosti populacije. 2.3 Določitev načrta vzorčenja Zanima nas povprečna teža bolnikov (µ) s hipertenzijo v starostni skupini 60 do 80 let. Težo bi radi ocenili na podlagi vzorca, jasno je, da bo teža precej različna pri moških (µ1) kot pri ženskah (µ2). Čas in denar, ki ju imamo na voljo za raziskavo, nam dopuščata vzorec velikosti 100. Vemo, da se v populaciji deleža moških in žensk s hipertenzijo razlikujeta, delež moških označimo z d. Zanima nas, kolikšen delež moških in kolikšen delež žensk naj naberemo v vzorec, da bo standardna napaka naše ocene najmanjša možna. Pri tem predpostavimo, da je standardni odklon teže moških k-krat večji od standardnega odklona teže žensk. 48 Verjetnost in statistika 2.3. Določitev načrta vzorčenja a) Zapišite nepristransko cenilko za populacijsko povprečje. Populacijsko povprečje µ lahko zapišemo kot µ = dµ1 + (1 − d)µ2. Ker velja E( ¯ X1) = µ1 in E( ¯ X2) = µ2, je nepristranska cenilka enaka M = d ¯ X1 + (1 − d) ¯ X2 b) Standardno napako ocene izrazite z velikostima podvzorcev (n1 je število moških, n2 število žensk). Ker je povprečje pri moških neodvisno od povprečja pri ženskah, velja σ2 σ2 var(M ) = d2var( ¯ X 1 2 1) + (1 − d)2var( ¯ X2) = d2 + (1 − d)2 n1 n2 c) Naj velja σ1 = kσ2. Pri kakšni razdelitvi vzorca je standardna napaka najmanjša? Izračunajte n1 za primera k = 1 in k = 2, vzemite, da je delež moških enak 0,7. Varianco vzorčnega povprečja zapišemo kot d2k2 (1 − d)2 var(M ) = σ2 + n1 n − n1 Minimizirati moramo torej izraz v oklepaju. Odvajamo po n1 in izenačimo z 0 d2k2 (1 − d)2 − + = 0 n2 (n − n 1 1)2 −(n − n1)2d2k2 + n2(1 − d)2 = 0 1 (d2k2 − (1 − d)2)n2 − 2nd2k2n 1 1 + n2d2k2 = 0 Rešimo kvadratno enačbo in dobimo rešitev (druga rešitev je večja od n in tako nesmiselna): ndk n1 = dk + 1 − d Za k = 1 dobimo rešitev n1 = nd, v konkretnem primeru torej n1 = 70 in n2 = 30. Če sta deleža moških in žensk na vzorcu enaka deležema v 49 Verjetnost in statistika 2.3. Določitev načrta vzorčenja populaciji, je standardna napaka najmanjša možna. Če je standardni odklon teže moških 2x večji od standardnega odklona žensk, dobimo n1 = 82,4, vzeti moramo torej ustrezno več posameznikov iz bolj variabilne skupine. d) Ali porazdelitev cenilke M konvergira proti normalni? V splošnem ne, to smo že pokazali v nalogi 1.1e). Če sta povprečji stra- tumov različni, bo porazdelitev cenilke odvisna od deleža d. Predlogi za vaje v R • Izmislite si smiselne vrednosti za parametre v nalogi ter generirajte po- datke. Grafično prikažite, kako se pri različnih izbirah velikosti vzorcev spreminja kvaliteta vaše ocene. • Poglejte si kaj se dogaja s porazdelitvijo cenilke za povprečje, ko se vzorec veča. Povzetek • Velikost vzorca ne vpliva na pristranskost ocene, je pa ključnega po- mena pri njeni standardni napaki. Pri stratificiranem vzorčenju na standardno napako vpliva velikost vzorcev in variabilnost v posameznih stratumih. Če je variabilnost v stratumih enaka, si najmanjšo standar- dno napako zagotovimo s proporcionalnim vzorčenjem, torej takim, pri katerem je delež enot izbranih v podvzorec enak deležu enot podsku- pine v populaciji. Če variabilnost v stratumih ni enaka, moramo iz bolj variabilnih pod- skupin vzeti ustrezno večji vzorec. • Če je stratumov več, je izpeljava nekoliko težja - minimizirati je po- trebno po več spremenljivkah hkrati. Izpeljava z vezanim ekstremom (omejitev n1 + . . . + nk = n) je podana v knjigi Rice (2009), razdelek 7.5.3. • Opisani vzorčni načrt je primer stratificiranega vzorčenja: populacija je razdeljena na več skupin (stratumov), vsak vzorec je sestavljen tako, da je vanj izbranih nekaj elementov vsakega stratuma. Vzorčna shema, pri kateri najprej naključno izberemo nekaj skupin in nato vzorčimo le iz njih, je primer večstopenjskega vzorčenja, srečali jo bomo v nalogi 2.6. 50 Verjetnost in statistika 2.4. Ocena kovariance 2.4 Ocena kovariance V nekem podjetju velikosti N so izvedli izobraževanje za naključen vzorec n zaposlenih. Ob koncu izobraževanja so novo znanje preverili s testom. Podjetje se želi odločiti, ali je smiselno uvesti izobraževanje za vse zaposlene, zato jih zanima povezanost med starostjo zaposlenega (Xi) in rezultatom na testu (Yi). Za vsakega posameznika iz vzorca imamo torej par slučajnih spremenljivk (Xi, Yi), i = 1 . . . n. a) Utemeljite, da je količina cov(Xi, Yj) za poljubna i = j enaka. Vzorčenje si lahko predstavljamo tako, da smo populacijo naključno ure- dili, nato pa v vzorec zajeli prvih n posameznikov. Ker imajo vsi vrstni redi enako verjetnost, bo na i-tem mestu z enako verjetnostjo katerikoli posameznik. Vsi pari (Xi, Yj) imajo tako enako porazdelitev in zato je enaka tudi kovarianca Xi in Yj. b) Naj bo γ = cov(Xi, Yi). Izračunajte kovarianco cov(Xi, Yj) za i = j. Uporabimo isto idejo kot pri prejšnji nalogi - vsota vseh vrednosti iz populacije je konstanta, zato velja N cov(Xi, Yj) = cov(Xi, Yi) + (N − 1)cov(Xi, Yj) = 0 j=1 Velja torej (za i = j) γ cov(Xi, Yj) = − N − 1 c) Kovarianca med spremenljivkama X in Y je definirana kot N N 1 1 cov(X, Y ) = cov(X1, Y1) = [(xi − µ)(yi − ν)] = xiyi − µν N N i=1 i=1 kjer smo z µ in ν označili povprečji µ = 1 N x N y N i=1 i in ν = 1 N i=1 i. 51 Verjetnost in statistika 2.4. Ocena kovariance S pomočjo vzorca bi kovarianco radi ocenili s cenilko n γ = c XiYi − n ¯ X ¯ Y i=1 Določite vrednost konstante c, da bo cenilka nepristranska. Pričakovana vrednost cenilke je n E(γ) = c E(XiYi) − nE( ¯ X ¯ Y ) (2.1) i=1 Zaradi simetrije je E(XiYi) = E(XjYj) za poljubna i in j. Vemo, da velja cov(Xi, Yi) = E(XiYi) − E(Xi)E(Yi) = E(XiYi) − µν Torej je E(XiYi) = µν + γ. Oglejmo si še drugi člen na desni strani (2.1): n n 1 1 E( ¯ X ¯ Y ) = E Xi Yj n n i=1 j=1 n n 1 = E XiYi + Xi Yj n2 i=1 j=1,i=j n n 1 = E(XiYi) + E(XiYj) n2 i=1 j=1,i=j n 1 = [E(XiYi) + (n − 1)E(XiYj)] n2 i=1 Uporabimo rezultat cov(Xi, Yj) = E(XiYj) − E(Xi)E(Yi) γ E(XiYj) = µν − N − 1 52 Verjetnost in statistika 2.4. Ocena kovariance in zato 1 −γ E( ¯ X ¯ Y ) = n µν + γ + (n − 1)(µν + ) n2 N − 1 1 (n − 1) = nµν + γ(1 − ) n N − 1 1 N − n = nµν + γ n N − 1 To vstavimo v enačbo (2.1) n 1 N − n E(γ) = c (µν + γ) − n nµν + γ n N − 1 i=1 N − n = c nµν + nγ − nµν − γ N − 1 N − n = c nγ − γ N − 1 nN − n N − n = cγ − N − 1 N − 1 N (n − 1) = cγ N − 1 c mora biti torej enak 1 N −1 . n−1 N d) Kako bi ocenili korelacijo? Kaj vemo o nepristranskosti te cenilke? Uporabimo izpeljane formule za oceno kovariance in varianc: cov(X, Y ) ρ = σXσY 1 N −1 n (Xi − ¯ X)(Yi − ¯ Y ) = n−1 N i=1 1 N −1 n (X N −1 n (Y n−1 N i=1 i − ¯ X)2 1 n−1 N i=1 i − ¯ Y )2 n (Xi − ¯ X)(Yi − ¯ Y ) = i=1 n (X (Y i=1 i − ¯ X)2 n i=1 i − ¯ Y )2 Nepristranskosti te ocene nismo dokazali, saj pričakovana vrednost kvo- cienta ni enaka kvocientu pričakovanih vrednosti. 53 Verjetnost in statistika 2.4. Ocena kovariance Predlogi za vaje v R • Ker ne vemo, ali je ocena korelacije pristranska ali ne, pristranskost preverimo s simulacijo. Vzamemo populacijo velikosti N = 300, vzorci naj bodo velikosti n = 10. Naj bo X starost porazdeljena enakomerno med 25 in 65, uspeh na testu pa negativno povezan s starostjo, tako da je v povprečju enak 100−starost (predpostavimo, da so odstopanja od tega povprečja razpršena s standardnim odklonom 20 in normalno porazdeljena). > set.seed(1) > xi <- runif(300)*40+25 #300 oseb starosti 25-65 let > yi <- 100 - xi + rnorm(300)*20 #rezultat testa za populacijo > cov(xi,yi) #kovarianca v populaciji [1] -136.8110 > cor(xi,yi) #korelacija v populaciji [1] -0.5207052 > runs <- 10000 #stevilo korakov simulacije > cova <- cora <- rep(NA,runs) #sem bomo zapisali rezultate > for(it in 1:runs){ #simulacija po korakih + inx <- sample(1:length(xi),size=10,replace=F) #vzorec n=10 + xa <- xi[inx] #pogledamo njihove starosti + ya <- yi[inx] #pogledamo njihove rezultate + cova[it] <- 1/9*299/300* + sum( (xa-mean(xa))*(ya-mean(ya))) #izracunamo kovarianco + cora[it] <- sum( (xa-mean(xa))*(ya-mean(ya)))/ + sqrt(sum( (xa-mean(xa))^2)*sum((ya-mean(ya))^2)) #korelacija + } > mean(cova) #povprecna kovarianca [1] -135.4745 > mean(cora) #povprecna korelacija [1] -0.5034081 • Vidimo, da sta obe vrednosti po absolutni vrednosti nekoliko manjši od populacijskih, preverimo, ali je odstopanje veliko glede na standardno napako, ki jo lahko pričakujemo pri takem številu simulacij. Zanima nas, ali povprečna kovarianca (mean(cova)) bistveno odstopa od prave vrednosti (cov(xi,yi)). Povprečna kovarianca je slučajna spremenljivka, če bomo vnovič pognali simulacijo (vseh 10000 kora- kov), bomo dobili drugo vrednost. Predpostavimo, da je približno nor- malno porazdeljena, ocenimo njeno varianco (varianca povprečja n i.i.d 54 Verjetnost in statistika 2.5. Enostavni vzorec, končna populacija II spremenljivk je varianca spremenljivk deljeno z n, pri nas je n število korakov simulacije). Ničelna domneva, ki jo preverjamo, je: H0: pov- prečna kovarianca je enaka populacijski vrednosti. Odstopanje od te ničelne domneve preverjamo s testom t. > (mean(cova)-cov(xi,yi))/sqrt(var(cova)/runs) [1] 1.509540 Ta rezultat je v okviru pričakovanj, saj smo teoretično pokazali, da je ocena kovariance nepristranska. Enako ponovimo za korelacijo: > (mean(cora)-cor(xi,yi))/sqrt(var(cora)/runs) [1] 6.66459 Odstopanje pri korelaciji je bistveno večje, verjamemo, da se v naši simulaciji ni zgodilo po naključju, temveč je ocena dejansko pristranska. Povzetek • Cilj naloge je bil oceniti povezanost med vrednostima dveh slučajnih spremenljivk na isti enoti neke populacije. Ker je populacija v našem primeru končna, je potrebno pri tej oceni upoštevati tudi kovarianco med vrednostmi na različnih enotah. Pokazali smo nepristranskost ocene kovariance, ne pa tudi nepristranskosti ocene korelacije. Za sled- njo se izkaže, da je pristranska, in sicer nekoliko podceni pravo vrednost (je bližje 0). 2.5 Enostavni vzorec iz končne populacije, še enkrat Vzemimo še enkrat enostavni slučajni vzorec velikosti n iz populacije N , vrednosti v populaciji označimo z xi; i = 1, . . . , N , populacijsko vrednost povprečja označimo z µ, variance pa z σ2. Definirajmo slučajno spremenljivko Ii = I[i je izbran v vzorec] in zapišimo cenilko populacijskega povprečja µ kot C = 1 N I n i=1 ixi. 55 Verjetnost in statistika 2.5. Enostavni vzorec, končna populacija II a) Koliko je vsota N I i=1 i? Kolikšna je verjetnost P (Ii = 1)? Vsota N I i=1 i = n, saj smo vzeli vzorec velikosti n. Izračunajmo še verjetnost, da bo izbran element i: izbiramo vzorce velikosti n in iz populacije velikosti N . Vseh možnih kombinacij je N , kakšno je število tistih vzorcev, v katerih je element n i? Pri teh vzorcih en element že poznamo, izmed ostalih N − 1 smo jih izbrali n − 1. Torej je iskana verjetnost enaka: N −1 n P (I n−1 i = 1) = = N N n b) Pokažite, da je cenilka nepristranska. Radi bi ocenili µ = 1 N x N i=1 i N 1 E(C) = E(Ii)xi n i=1 Ker lahko Ii zavzame le vrednosti 0 in 1, je E(Ii) = P (Ii = 1) = n N (vzorec je slučajen, zato so verjetnosti za vse i enake), zato dobimo N N 1 n 1 E(C) = xi = xi = µ n N N i=1 i=1 c) Izračunajte var(Ii) in cov(Ii, Ij). Spremenljivka Ii je Bernoullijeva z verjetnostjo P (Ii = 1) = n . Njena N varianca je zato enaka n n n N − n var(Ii) = (1 − ) = N N N N Kovarianco izračunamo tako, da upoštevamo, da je cov(I1, I1 + . . . IN ) = cov(I1, n) = 0 in cov(Ii, Ij) je enaka za vsak i = j: n N −n n(N − n) cov(I N N i, Ij ) = − = − N − 1 N 2(N − 1) 56 Verjetnost in statistika 2.5. Enostavni vzorec, končna populacija II d) Pokažite še, da je varianca tako zapisane cenilke enaka var(C) = σ2 N−n . n N −1 N N 1 var(C) = cov( Iixi, Ijxj) n2 i=1 j=1 N N 1 = cov(Iixi, Ijxj) n2 i=1 j=1 N N 1 = cov(Iixi, Iixi) + cov(xiIi, Ijxj) n2 i=1 j=1,j=i N N 1 = x2cov(Ii, Ii) + xixjcov(Ii, Ij) n2 i i=1 j=1,j=i Populacijska varianca je definirana kot: N N 1 1 σ2 = (xi − µ)2 = x2 − µ2 N N i i=1 i=1 zato velja N x2 = N (σ2 + µ2). V izpeljavi variance bomo potrebovali i=1 i še naslednji rezultat: N N N N N N N ( xi)2 = x2 + x x x2 i ixj ⇒ ixj = N 2µ2 − i i=1 i=1 i=1 j=1,j=i i=1 j=1,j=i i=1 57 Verjetnost in statistika 2.6. Večstopenjsko vzorčenje Sedaj to uporabimo in dobimo N N 1 var(I var(C) = x2var(I i) i) − xixj n2 i N − 1 i=1 j=1,j=i N N var(I = i) (N − 1) x2 − (N 2µ2 − x2) n2(N − 1) i i i=1 i=1 N var(I = i) N x2 − N 2µ2 n2(N − 1) i i=1 N N 2var(I 1 = i) x2 − µ2 n2(N − 1) N i i=1 N 2var(I = i) σ2 n2(N − 1) N 2σ2 n N − n = n2(N − 1) N N σ2 N − n = n N − 1 Povzetek • V tej nalogi smo ponovili rezultate iz naloge 2.2, vendar z nekoliko drugačnim zapisom vzorčnega povprečja. Medtem ko so bile izpeljave na tem primeru zaradi tega nekoliko zapletenejše, se bo ta pristop iz- kazal za ključnega v naslednji nalogi, ko bo populacija razdeljena na več skupin. Ta naloga tako služi predvsem kot nekoliko bolj pregleden uvod v naslednjo. 2.6 Večstopenjsko vzorčenje Oceniti želimo dosežek ljubljanskih sedmošolcev na nekem testu znanja, ki ga izvajajo v več državah. Populacijo N = 2800 učencev te starosti bomo vzorčili po šolah (K = 46). V vzorec bomo najprej slučajno (in neodvisno od števila Ni sedmošolcev na šoli i) vzorčili k = 10 šol, nato pa bomo na vsaki šoli izbrali vzorec n = 15 učencev. Naj µ označuje populacijsko povprečje dosežka na testu, µi pa naj bo povprečje za vsako šolo posebej. Vzorčenje znotraj šol je neodvisno od vzorčenja na prvem koraku. 58 Verjetnost in statistika 2.6. Večstopenjsko vzorčenje a) Zapišite nepristransko cenilko za µ. Najprej izrazimo µ s povprečji šol, torej µi. Naj bo xij vrednost j-tega učenca na i-ti šoli. Velja K N K 1 i 1 µ = xij = Ni · µi (2.2) N N i=1 j=1 i=1 Označimo ocenjeno povprečje vsake šole z ¯ Xi, Ii pa naj bo indikatorska spremenljivka, ki je enaka 1, če je šola izbrana v vzorec. Naša cenilka naj bo enaka K ¯ X = c ¯ iIiXi i=1 Določiti moramo vrednost konstante ci, tako da bo cenilka nepristranska. Upoštevamo, da smo na vsaki šoli vzeli naključni vzorec, in zato velja E( ¯ Xi) = µi. Ker je vzorčenje na drugem koraku neodvisno od vzorčenja na prvem, velja E(I ¯ iXi) = E(Ii)E( ¯ Xi). Ker smo na prvem koraku vzorčili vse šole z enako verjetnostjo, je E(Ii) = k za vsak i. Uporabimo vse K našteto in dobimo K K E( ¯ X) = c ¯ iE(IiXi) = ciE(Ii)E( ¯ Xi) i=1 i=1 K k = ci µi K i=1 Zaradi (2.2) mora veljati c k i = Ni , zato je naša cenilka enaka K N K ¯ K 1 X = N ¯ iIiXi N k i=1 b) Kako bi ocenili populacijsko povprečje, če bi imele vse šole enako število učencev L? Ker velja N = K N i=1 i, za enake Ni = L velja N = K L, in zato K K ¯ 1 1 1 X = LI ¯ ¯ iXi = IiXi L k k i=1 i=1 59 Verjetnost in statistika 2.6. Večstopenjsko vzorčenje c) Ali je za nepristranskost pomembno, koliko učencev z vsake šole vzamete? Ne, ¯ Xi je nepristranska cenilka µi ne glede na velikost vzorca. Seveda pa velikost vzorca vpliva na standardno napako te cenilke. d) Zapišite varianco cenilke s pomočjo varianc in kovarianc K K 1 var( ¯ X) = var( N ¯ iIiXi) N k i=1 K K 2 K = N 2var(I ¯ ¯ ¯ iXi) + NiNjcov(IiXi, IjXj) N k i i=1 j=1,i=j e) Označimo varianco znotraj vsake šole z σ2 = 1 Ni (x wi N ij − µi)2. Kaj je i j=1 var(I ¯ ¯ ¯ iXi) in kaj cov(IiXi, Ij Xj )? Uporabimo, da je vzorčenje na drugem koraku neodvisno od vzorčenja na prvem in da je I2 = I i i (12 = 1, 02 = 0): var(I ¯ ¯ ¯ iXi) = E(I2X2) − E(I X ) − E(I i i i i)2 = E(Ii)E( ¯ X2i i)2E( ¯ Xi)2 k k 2 = E( ¯ X2) − µ2 K i K i Upoštevamo še, da je E( ¯ X2) = var( ¯ X Ni−n + µ2, in i i) + E( ¯ Xi)2 = σ2wi n N i i−1 dobimo k σ2 N k 2 var(I ¯ wi i − n iXi) = + µ2 − µ2 K n N i i i − 1 K k(K − k) k σ2 N = µ2 + wi i − n i K2 K n Ni − 1 Sedaj izrazimo še kovarianco: cov(I ¯ ¯ ¯ ¯ ¯ ¯ iXi, Ij Xj ) = E(IiIj XiXj ) − E(IiXi)E(Ij Xj ) 60 Verjetnost in statistika 2.6. Večstopenjsko vzorčenje Upoštevamo neodvisnost vzorčenja na prvem in drugem koraku in dejstvo, da je povprečje na eni šoli neodvisno od povprečja druge šole: cov(I ¯ ¯ iXi, Ij Xj ) = E(IiIj)µiµj − E(Ii)E(Ij)µiµj k(K − k) = µiµjcov(Ii, Ij) = −µiµj K2(K − 1) f) Izpeljite formulo za varianco cenilke v primeru, ko so vse vrednosti Ni enake L in je varianca znotraj šole enaka za vse šole, varianco med šolami označite z σ2. b K 1 2 K var( ¯ X) = L2var(I ¯ ¯ ¯ iXi) + L2cov(IiXi, IjXj) Lk i=1 i=1,i=j 1 2 K k(K − k) k σ2 L − n = µ2 + w k i K2 K n L − 1 i=1 K k(K − k) − µiµj K2(K − 1) j=1,i=j Uporabimo izraz za varianco med šolami K K K 1 1 1 σ2 = [µ µ2 − 2µµ µ2 − µ2 b i − µ]2 = i − µ2 = K K i K i i=1 i=1 i=1 in podobno kot pri nalogi 2.5 K K K µiµj = µ2K2 − µ2i i=1 j=1,i=j i=1 61 Verjetnost in statistika 2.6. Večstopenjsko vzorčenje K K K k(K − k) k(K − k) µ2 − µ i iµj K2 K2(K − 1) i=1 i=1 j=1,i=j K K K k(K − k) = (K − 1) µ2 − µiµj K2(K − 1) i i=1 i=1 j=1,i=j K K k(K − k) = (K − 1) µ2 − (µ2K2 − µ2) K2(K − 1) i i i=1 i=1 K k(K − k) = K µ2 − µ2K2 K2(K − 1) i i=1 K k(K − k) 1 = µ2 − µ2 K − 1 K i i=1 k(K − k) = σ2 K − 1 b in zato K 1 k(K − k) k σ2 L − n var( ¯ X) = σ2 + w k2 (K − 1) b K n L − 1 i=1 1 k(K − k) K k σ2 L − n = σ2 + w k2 (K − 1) b k2 K n L − 1 σ2 K − k 1 σ2 L − n = b + w k (K − 1) k n L − 1 g) Ali lahko porazdelitev cenilke za povprečje aproksimiramo z normalno porazdelitvijo? Da bi porazdelitev cenilke konvergirali proti neskončnosti ni dovolj le naraščanje velikosti vzorca - nujen pogoj je, da gre proti neskončnosti tudi število izbranih skupin. Predlogi za vaje v R • Preverite rezultate naloge z R. 62 Verjetnost in statistika 2.6. Večstopenjsko vzorčenje • S simulacijami si oglejte, kaj se dogaja s porazdelitvijo cenilke pov- prečja, ko gre velikost vzorca proti neskončnosti. Povzetek • Ko je populacija razdeljena na skupine, ki imajo različna povprečja, povprečje enot, zajetih v vzorec, ni več nujno nepristranska cenilka populacijskega povprečja - v ta namen je povprečja posameznih skupin potrebno ustrezno utežiti. • Možne so različne sheme vzorčenja, v katerih imajo lahko nepristranske cenilke pri enakem številu enot različno standardno napako. Pogosto bo tako cilj poiskati shemo vzorčenja, ki nam pri enaki skupni veli- kosti vzorca da najbolj natančen rezultat, torej najmanjšo standardno napako. • Oglejmo si še enkrat izraz za varianco cenilke in ga skušajmo interpre- tirati. Denimo, da imamo veliko število šol in veliko število učencev na vsaki šoli (popravka za končno populacijo postaneta zanemarljiva), in upoštevajmo, da so vse šole enako velike (enako prispevajo k skupnemu povprečju) in da je varianca znotraj vseh šol enaka. Dobljeni izraz za varianco cenilke poenostavi v σ2 σ2 var( ¯ X) = b + w k kn Vidimo, da za varianco cenilke pomembnejšo vlogo igra varianca med povprečji šol, medtem ko del, ki izhaja iz variabilnosti znotraj šol, hitro postane zelo majhen. Pri dani velikosti vzorca bo cenilka torej imela najmanjšo varianco ta- krat, ko bomo vzeli le po enega učenca z vsake šole. V tem primeru je vzorec sestavljen iz neodvisnih slučajnih spremenljivk in tako se noben del informacije ne podvaja zaradi odvisnosti. Je pa iz takega vzorca se- veda nemogoče ocenjevati varianco znotraj skupin, pa tudi o povprečjih skupin vemo zelo malo. 63 Poglavje 3 Ocenjevanje parametrov - metoda največjega verjetja V prejšnjem poglavju smo intuitivno smiselno cenilko za povprečje ali va- rianco določili kar neposredno, tako da smo sledili definiciji teh količin v populaciji. V splošnem bo smiselno cenilko pogosto precej težje določiti. Še pomembnejše kot le imeti ustrezno cenilko je tudi poznavanje njene poraz- delitve. V prejšnjem poglavju smo za aproksimacijo porazdelitve uporabili centralni limitni izrek - cenilke, ki so nas zanimale, se je namreč dalo izraziti kot vsote neodvisnih enako porazdeljenih slučajnih spremenljivk, kar pa v splošnem seveda ne bo vedno res. Metoda največjega verjetja (Rice, 2009, razdelek 8.5) je generična metoda, ki pomaga izpeljati smiselno cenilko, katere asimptotska porazdelitev je znana. Z večanjem velikosti vzorca n proti neskončnosti se porazdelitev cenilke po metodi največjega verjetja približuje normalni porazdelitvi. Vrednost te ce- nilke z večanjem vzorca konvergira k populacijski vrednosti, ki jo ocenjujemo, cenilka po metodi največjega verjetja je torej dosledna (vendar na majhnih vzorcih ne nujno nepristranska). Hkrati je z metodo podan tudi postopek izpeljave variance, s čimer je asimptotska porazdelitev določena. Ključno vlogo v teoritični podlagi te metode igrata prvi in drugi odvod lo- garitma verjetja (Rice, 2009, razdelek 8.5.2). Prvi odvod imenujemo funk- cija zbira (‘score function’), z njegovo pomočjo izpeljemo cenilko, negativno pričakovano vrednost drugega odvoda imenujemo Fisherjeva informacija in je del formule za varianco. Poglavje pričnemo z iskanjem cenilke za populacijski delež. Na tem primeru pokažemo, da metoda največjega verjetja predlaga enostavno cenilko, do ka- 64 Verjetnost in statistika 3.1. Ocenjevanje deleža tere bi tudi sicer prišli po intuitivni poti - na ta način na primeru skušamo razumeti smiselnost metode. Druga naloga v tem poglavju obravnava ravno tako zelo pogosto uporabljan problem - po metodi največjega verjetja izpelje cenilki za parametra v linearni regresiji (Rice, 2009, poglavje 14). Ta naloga tako služi za uvod v naslednja poglavja kot tudi za primer, v katerem hkrati ocenjujemo več kot en parameter. 3.1 Ocenjevanje deleža Naj bodo x1, . . . xn neodvisne realizacije Bernoullijevo porazdeljene slučajne spremenljivke X. Radi bi ocenili parameter p. a) Recimo, da je n = 5 in da smo dobili naslednjih 5 vrednosti: 1,0,1,1,1. Kolikšna bi bila verjetnost tega dogodka, če bi bil p = 0,2? Kaj pa za p = 0,75? Narišite krivuljo verjetnosti tega dogodka glede na p. Kako bi izračunali njen vrh? Verjetnost dogodka izračunamo kot 0,240,81, torej pk(1 − p)n−k, kjer je k število enic. Označimo z A dogodek A = {X1 = 1, X2 = 0, X3 = 1, X4 = 1, X5 = 1}. Za p = 0,2 dobimo P (A) = 0,00128, za p = 0,75 dobimo P (A) = 0,079. Narišemo krivuljo za vrednosti p med 0 in 1: 0.08 0.06 y 0.04 0.02 0.00 0.0 0.2 0.4 0.6 0.8 1.0 p Slika 3.1: Verjetnost opaženega dogodka glede na p. Vrh funkcije lahko poiščemo z odvajanjem - odvajamo funkcijo pk(1−p)n−k po p in izenačimo z 0 (lokalni maksimum). Vrh ni odvisen od vrstnih 65 Verjetnost in statistika 3.1. Ocenjevanje deleža redov. V našem primeru je vrh funkcije dosežen pri p = 4 . 5 b) Podatke, ki jih dobimo na nekem vzorcu, označimo z x1, . . . , xn (v zgor- njem primeru je bil n = 5, x1 = 1 in x2 = 0). Za vsako enoto zapišite P (Xi = xi|p), torej verjetnost, da se je zgodil dogodek, ki smo ga videli. Zapišite funkcijo verjetja. P (Xi = xi|p) = pxi(1 − p)1−xi Funkcija verjetja je produkt posameznih verjetnosti (predpostavili smo, da so slučajne spremenljivke Xi neodvisne), torej n L(p, x) = P (X1 = x1, . . . , Xn = xn|p) = pxi(1 − p)1−xi i=1 n = p x x i=1 i (1 − p)n− ni=1 i c) Poiščite oceno za p po metodi največjega verjetja. Ker je logaritem monotona funkcija, lahko namesto lokalnega maksimuma te funkcije gledamo raje maksimum logaritma: n n log L(p, x) = xi log(p) + (n − xi) log(1 − p) i=1 i=1 ∂ log L(p, x) n xi n − n xi = i=1 − i=1 ∂p p 1 − p n xi − p n xi − p(n − n xi) = i=1 i=1 i=1 p(1 − p) n xi − pn = i=1 p(1 − p) Odvod logaritma verjetja bo enak 0 pri ˆ pn = n x i=1 i. Ocena po metodi največjega verjetja je torej ˆ p = 1 n x n i=1 i. Ocena je ravno delež enic v vzorcu. 66 Verjetnost in statistika 3.1. Ocenjevanje deleža d) Ali je ocena nepristranska? Metoda največjega verjetja zagotavlja le doslednost (nepristranost, ko gre n → ∞), v našemo primeru dobimo n n n 1 1 1 E(ˆ p) = E( xi) = E(xi) = p = p n n n i=1 i=1 i=1 V našem primeru je torej ocena nepristranska. e) Zapišite oceno standardne napake Varianca ocene je enaka 1 I(p)−1, kjer I(p) označuje Fisherjevo informa- n cijo: ∂2 ∂ 2 I(p) = −E log(f (X, p)) = E log(f (X, p)) ∂p2 ∂p V našem primeru sta izračuna po obeh formulah enako težka, uporabimo prvo formulo: f (X|p) = pX (1 − p)1−X ∂2 I(p) = −E log(f (X|p)) ∂p2 ∂2 = −E (X log p + (1 − X) log(1 − p)) ∂p2 ∂ X 1 − X = −E − ∂p p 1 − p ∂ (1 − p)X − (1 − X)p = −E ∂p p(1 − p) ∂ X − p = −E ∂p p(1 − p) p(1 − p)(−1) − (1 − 2p)(X − p) = −E p2(1 − p)2 −p + p2 − X + 2pX + p − 2p2 = −E p2(1 − p)2 −p2 − X + 2pX = −E p2(1 − p)2 67 Verjetnost in statistika 3.1. Ocenjevanje deleža Pri računanju pričakovane vrednosti upoštevamo, da je E(X) = p, ker je X le v števcu, dobimo −p2 − X + 2pX I(p) = −E p2(1 − p)2 −p + p2 = − p2(1 − p)2 1 = p(1 − p) f) Oceniti želimo delež volilcev nekega kandidata. Na vzorcu n = 500 zanj glasuje 29 % volilcev. Podajte 95% interval zaupanja za to oceno. Vzorčna ocena je ˆ p = 0,29. Standardno napako (torej standardni od- klon cenilke) na vzorcu ocenimo s pomočjo ˆ p, ocena standardne napake je torej enaka 1 ˆ p(1 − ˆ p) SE = = = 0,02 nI(ˆ p) n Teorija nam pove, da je lahko porazdelitev kvocienta p−ˆp aproksimiramo SE z normalno porazdelitvijo, 95% interval zaupanja je enak [0,25, 0,33]. Predlogi za vaje v R • Z R narišite sliko 3.1: > p <- seq(0,1,length=100) #za 100 vrednosi p med 0 in 1 > y <- p^4*(1-p) #za vsako vrednost izracunam verjetnost > plot(p,y,type="l") #narisem in povezem s krivuljo • Generirajte vzorec velikosti 500, v katerem ima vsak posameznik verje- tnost 0,3, da glasuje za nekega kandidata. Ocenite verjetnost z deležem na vzorcu. Ponovite poskus 1000x in si oglejte porazdelitev vzorčnih ocen. • Na vsakem vzorcu ocenjenemu deležu dodajte še 95% interval zaupanja. Kolikšen je delež vzorcev, pri katerih interval zaupanja zajema pravo vrednost (0,3)? Temu deležu pravimo pokritje (angl. ’coverage’). Pozor - ta delež ni nujno enak 95%, saj nam teorija da le približek porazdelitvi, na manjših vzorcih bo zato ta odstotek lahko precej odstopal. 68 Verjetnost in statistika 3.2. Povezanost dveh spremenljivk Povzetek • V nalogi smo za preprost primer narisali funkcijo verjetja in intuitivno utemeljili, zakaj je točka, v kateri ta funkcija doseže vrh, smiselna ocena parametra. • Ocena deleža v populaciji po metodi največjega verjetja je kar delež na vzorcu. Pokazali smo, da je ta ocena nepristranska. • Varianca ocene deleža je enaka ˆ p(1− ˆ p) . Odvisna je od velikosti vzorca n in od dejanske vrednosti deleža, ki ga ocenjujemo - napaka bo največja pri ocenjevanju deležev blizu polovice. 3.2 Povezanost dveh spremenljivk Zanima nas, kako je prihodek podjetja v neki panogi odvisen od števila za- poslenih. Predpostavimo, da je prihodek podjetja normalno porazdeljen s povprečjem β0 + β1X, kjer je X logaritem števila zaposlenih. Denimo, da imamo podatke o številu zaposlenih in prihodku za vzorec podjetij, radi bi ocenili parametra β0 in β1. a) Zapišite gostoto porazdelitve prihodka podjetja, če vemo, da je varianca enaka σ2. Predpostavljamo, da je Y ∼ N (β0 + β1X, σ2), torej 1 f (Y, X|β0, β1, σ) = √ e− (Y −β0−β1X)2 2σ2 2πσ b) Zapišite funkcijo verjetja. Kaj je funkcija, ki jo moramo maksimizirati? Dani so podatki (xi, yi), i = 1, . . . , n. n 1 L(y, x, β0, β1, σ) = √ e− (yi−β0−β1xi)2 2σ2 2πσ i=1 1 n i=1(yi−β0−β1xi)2 = √ e− 2σ2 ( 2πσ)n 69 Verjetnost in statistika 3.2. Povezanost dveh spremenljivk Logaritem te funkcije je n n (yi − β0 − β1xi)2 log L(y, x, β i=1 0, β1, σ) = − log(2πσ2) − 2 2σ2 Ker nas zanimata le parametra β0 in β1, je prvi del funkcije konstanta, maksimizirati je potrebno le izraz n − (yi − β0 − β1xi)2 i=1 Izraz yi−(β0+β1xi) predstavlja razdaljo med točkama T (xi, yi) in T (xi, β0+ β1xi), to vrednost imenujemo ostanek (razdalja točke od premice). Oceni za parametra β0 in β1 določata premico, ki se najbolje prilega podatkom v smislu, da je vsota kvadriranih ostankov točk od premice najmanjša možna. To oceno zato imenujemo ocena po metodi najmanjših kvadratov (Rice, 2009, razdelek 14.1). c) Izračunajte oceni β0 in β1 po metodi največjega verjetja. Najprej za β0: n ∂ (yi − β0 − β1xi)2 ∂β0 i=1 n = −2 (yi − β0 − β1xi) i=1 Če zgornji izraz izenačimo z 0, dobimo (izraz je enak nič za posebni vre- dnosti β0 in β1, ki ju označimo s strešico) n n −2 yi − nβ0 − β1 xi = 0 i=1 i=1 n n 1 β0 = yi − β1 xi n i=1 i=1 70 Verjetnost in statistika 3.2. Povezanost dveh spremenljivk Sedaj odvajamo še po β1: n ∂ (yi − β0 − β1xi)2 ∂β1 i=1 n = −2 xi (yi − β0 − β1xi) i=1 n n n = −2 xiyi − β0 xi − β1 x2i i=1 i=1 i=1 Če zgornji izraz izenačimo z 0, dobimo n n xiyi − β0 xi β i=1 i=1 1 = n x2i i=1 Združimo obe izpeljavi in (po malce premetavanja členov) dobimo n n n n xiyi − xi yi β i=1 i=1 i=1 1 = n n n x2 − ( x i i)2 i=1 i=1 n n n n x2 y x x i i − i iyi β i=1 i=1 i=1 i=1 0 = n n n x2 − ( x i i)2 i=1 i=1 d) Izračunajte standardno napako za obe oceni. Za Fisherjevo matriko informacije moramo izračunati druge odvode. Lo- garitem funkcije verjetja je enak 1 (Y − β log f (Y, X|β 0 − β1X )2 0, β1, σ) = − log(2πσ2) − 2 2σ2 71 Verjetnost in statistika 3.2. Povezanost dveh spremenljivk Prva odvoda sta enaka ∂ 1 log f (Y, X|β0, β1, σ) = (Y − β0 − β1X) ∂β0 σ2 ∂ X log f (Y, X|β0, β1, σ) = (Y − β0 − β1X) ∂β1 σ2 Drugi odvodi so potem ∂2 1 log f (Y, X|β0, β1, σ) = − ∂β2 σ2 0 ∂2 X2 log f (Y, X|β0, β1, σ) = − ∂β2 σ2 1 ∂2 X log f (Y, X|β0, β1, σ) = − ∂β1β0 σ2 Členi Fisherjeve matrike informacije so negativne pričakovane vrednosti drugih odvodov. Ker pričakovane vrednosti X oziroma X2 ne poznamo, ju ocenimo iz podatkov: 1 ¯ x  1 I(β n 0, β1) =   σ2 ¯ x 1 x2 n i i=1 Inverz te matrike je potem  n  σ2 1 x2 −¯ x I−1(β n i 0, β1) = n  i=1  1 x2 − ¯ x2 −¯ x 1 n i i=1 in zato n σ2 1 x2 I−1 1 n i var(β 11 i=1 0) = = n n n 1 x2 − ¯ x2 n i i=1 n σ2 x2i = i=1 n n n x2 − ( x i i)2 i=1 i=1 72 Verjetnost in statistika 3.2. Povezanost dveh spremenljivk ter I−1 1 σ2 var(β 22 1) = = n n n 1 x2 − ¯ x2 n i i=1 nσ2 = n n n x2 − ( x i i)2 i=1 i=1 Povzetek • Če je pogojna porazdelitev spremenljivke Y glede na X normalna, je ocena parametrov v linearni regresiji po metodi največjega verjetja enaka oceni po metodi najmanjših kvadratov. • Metoda najmanjših kvadratov je intuitiven pristop, ki zagotavlja smi- selno premico. Da bi lahko uporabili to metodo, moramo definirati model in s tem ostanke, ni pa nam potrebno narediti nobenih predpo- stavk o porazdelitvi ostankov. Metoda nam da cenilke koeficientov, ne izvemo pa ničesar o njihovi standardni napaki oziroma porazdelitvi. V posameznih primerih bo navkljub temu porazdelitev dokaj preprosto najti, saj so cenilke koeficientov linearne kombinacije spremenljivke Y . Nasprotno metoda največjega verjetja zahteva predpostavke o poraz- delitvi ostankov ter nato tudi poda asimptotsko porazdelitev teh ocen glede na pravo vrednost parametrov v populaciji. 73 Poglavje 4 Preizkušanje domnev Statistično sklepanje, opisano v prejšnjih poglavjih je bilo osredotočeno, na ocenjevanje vrednosti v populaciji s pomočjo vzorca. V tem razdelku bomo s pomočjo vzorca preverjali resničnost trditev o populaciji. Frekventistična statistika pri tovrstnem preizkušanju domnev temelji na Ney- man-Pearsonovi paradigmi (Rice, 2009, razdelek 9.2), ki predstavlja način odločanja med domnevami. Pri tem osrednjo vlogo igra ena izmed domnev, imenujemo jo ničelna domneva, ki naj bi predstavljala naše dosedanje znan- je. Od nje smo se pripravljeni oddaljiti le, če je verjetnost, da vzorec izhaja iz populacije, v kateri ta domneva drži, zelo majhna. Da bi o tem lahko presojali, definiramo smiselno testno statistiko, ki predstavlja nek povzetek informacije o ničelni domnevi, ki nam jo ponuja vzorec. Če poznamo porazdelitev te testne statistike pod ničelno domnevo, lahko potem presojamo o ‘nenavadnosti’ njene vrednosti na našem vzorcu. Porazdelitev testne statistike pod ničelno domnevo nam bo omogočila postavitev mej, v katerih se z veliko verjetnostjo nahaja vrednost testne statistike za vzorce, ki izhajajo iz populacije, v kateri ničelna domneva velja. Če bodo meje na vzorcu presežene, bomo govorili o statistično značilnem rezultatu in ničelno domnevo zavrnili v korist alternativne. Pri tovrstnem odločanju lahko zagrešimo dve vrsti napak. Napaka prve vrste je verjetnost, da zavrnemo ničelno domnevo, četudi v populaciji drži. Velikost te napake določimo sami, glede na to kako veliko verjetnost te napake si bomo dovolili, bomo postavili meje oziroma določili območje zavrnitve. Največjo dovoljeno verjetnost napake bomo imenovali stopnja značilnosti in jo označili z α. Na- paka druge vrste (označili jo bomo z β) je verjetnost, da ničelne domneve na 74 Verjetnost in statistika 4.1. Preizkušanje domnev podlagi vzorca ne bomo zavrnili, četudi bi jo morali. O velikosti te napake lahko sklepamo, če poznamo porazdelitev testne statistike pod alternativno domnevo. Na primerih si bomo ogledali, kaj vse vpliva na njeno velikost. Pogosto bomo namesto o napaki druge vrste govorili o nasprotni verjetnosti, torej verjetnosti, da uspemo zavrniti ničelno domnevo, kadar vzorec dejansko ne izhaja iz populacije, v kateri bi ničelna domneva veljala. To verjetnost bomo imenovali moč testa. Poglavje pričnemo s preprostim primeom s pomočjo katerega uvedemo osnovne pojme in ideje statističnega preizkušanja domnev. Sledi konkre- ten primer iz prakse, ki nam pokaže, kako pomemben je lahko izračun moči pred začetkom raziskave. Nato se vrnemo na nekoliko splošnejši primer, pri katerem se na podlagi vzorca želimo odločati med dvema domnevama o po- pulaciji, a imamo vsaj za začetek domnevi za enakovredni. Odločanje poteka tako, da na vzorcu definiramo smiselno testno statistiko, katere vrednosti bodo ‘dokazi’ v prid eni ali drugi domnevi. Odločali se bomo glede na ver- jetnost, torej je nujno poznati porazdelitev te testne statistike, ki pa bo pod različnima domnevama seveda različna. Nalogo bomo nato še posplošili in si pogledali primer uporabe centralnega limitnega izreka pri iskanju porazdeli- tve testne statistike. Poglavje bomo zaključili s posplošenim testom razmerja verjetij (Rice, 2009, razdelek 9.4). Podobno kot je metoda največjega verjetja pomenila generičen pristop k iskanju smiselne cenilke, ki hkrati podaja tudi njeno asimptotsko porazdelitev, tudi posplošeni test razmerja verjetij pomeni generičen pristop, ki nam predlaga smiselno testno statistiko in poda njeno asimptotsko porazdelitev. 4.1 Osnovni pojmi pri statističnem preizku- šanju domnev Želimo preveriti, ali je kovanec pošten. Naredili smo poizkus, kjer smo 10- krat vrgli kovanec in pri tem 7-krat dobili grb. a) Zapišite ničelno domnevo za vaš primer. Ali je ničelna domneva eno- stavna ali sestavljena? Zapišite testno statistiko, označite jo z X - kakšna je njena porazdelitev pod ničelno domnevo? 75 Verjetnost in statistika 4.1. Preizkušanje domnev H0 : p = 0,5. S tem je porazdelitev slučajne spremenljivke pod ničelno domnevo natanko določena, zato pravimo, da je ničelna domneva eno- stavna. Testna statistika X je število grbov: porazdelitev testne statistike pod ničelno domnevo je binomska B(10, 0,5). b) Denimo, da je vaša alternativna domneva HA : p > 0,5. Ali je ta domneva enostavna ali sestavljena? Pri kakšnih vrednostih X boste zavrnili ničelno domnevo v prid alternativni? Ali je domneva enostranska ali dvostranska? Domneva je sestavljena, saj zajema več vrednosti istega parametra - tudi če vemo, da domneva drži, še ne poznamo porazdelitve. Zavračali bomo pri velikih vrednostih X. Domneva je enostranska, saj nas bo zanimal le desni rep porazdelitve. c) V našem primeru je X = 7. Kolikšna je verjetnost, da se na vzorcu zgodi ta dogodek, če ničelna domneva drži? Če ničelna domneva drži (p = 0,5), je P0(X = 7) = 0,117. d) Denimo, da je območje zavrnitve sestavljeno iz vrednosti {10}. Kolikšna je stopnja značilnosti α v tem primeru? Kolikšna je stopnja značilnosti, če je območje zavrnitve sestavljeno iz vrednosti {6,7,8,9,10}? Če ničelna domneva drži (p = 0,5), je P0(X = 10) = 0,001. Če ničelna domneva drži (p = 0,5), je P0(X ≥ 6) = 0,377. e) Določite območje zavrnitve pri stopnji značilnosti α = 0,05. Ali lahko na podlagi dobljenih podatkov zavrnete ničelno domnevo pri stopnji značilnosti α = 0,05? Največja vrednost k, tako da velja P0(X ≥ k) < 0,05, je enaka 9. Stop- nja značilnosti je v tem primeru enaka 0,01. Ničelne domneve seveda ne moremo zavrniti. f) Kolikšna je moč testa pri tej vrednosti α, če predpostavimo, da je prava vrednost parametra p = 0,6? Kaj pa pri p = 0,7? Kolikšna je v teh primerih napaka druge vrste? Moč testa je majhna - 0,046 oziroma 0,149. Pri tako majhnem vzorcu 76 Verjetnost in statistika 4.1. Preizkušanje domnev in majhni stopnji značilnosti bomo ničelno domnevo le redko uspešno za- vrnili. Napako druge vrste izračunamo kot 1 - moč. g) Predpostavite, da je vaša alternativna domneva enaka HA : p = 0,5. Ali je ta domneva enostavna ali sestavljena? Ali je domneva enostranska ali dvostranska? Alternativna domneva je še vedno sestavljena, zdaj je tudi dvostranska. h) Kakšno bo sedaj območje zavrnitve, če želite, da je α ≤ 0,05? Kolikšna natanko bo stopnja značilnosti za to območje? Območje zavrnitve bo sestavljeno iz vrednosti {0,1,9,10}. Stopnja značilnosti za to območje je enaka α = 0,02. i) Izračunajte še moč testa v tem primeru. Pri moči testa moramo upoštevati, da bomo sedaj ničelno domnevo zavr- nili tudi če bo X enak 0 ali 1. Ker pa je verjetnost teh dveh vrednosti za p = 0,6 oz. p = 0,7 zelo majhna, se moč praktično ne spremeni. Tabela 4.1: Kumulativne verjetnosti P (X ≤ k|p) za binomsko porazdelitev pri n = 10. p\k 0 1 2 3 4 5 6 7 8 9 10 0,5 0,001 0,011 0,055 0,172 0,377 0,623 0,828 0,945 0,989 0,999 1 0,6 0,000 0,002 0,012 0,055 0,166 0,367 0,618 0,833 0,954 0,994 1 0,7 0,000 0,000 0,002 0,011 0,047 0,150 0,350 0,617 0,851 0,972 1 Predlogi za vaje v R • 1000x ponovite poskus, v katerem po 10x mečete kovanec. Oglejte si verjetnost zavrnitve za posamezno območje. • Spremenite verjetnost, s katero pade grb, in si oglejte moč testa. • Povečajte vzorec in si oglejte, kako se spreminja moč testa. 77 Verjetnost in statistika 4.2. Moč testa Povzetek • V nalogi smo se z enostavnim primerom sprehodili čez osnovne pojme preizkušanja domnev, kot so: ničelna domneva, testna statistika, eno- stavna oziroma sestavljena domneva, stopnja značilnosti, zavrnitveno območje, enostranska oziroma dvostranska alternativna domneva, moč testa. • Neyman-Pearsonova lema (Rice, 2009, stran 332) pove, da ima raz- merje verjetij pri dani stopnji značilnosti α največjo moč izmed vseh testnih statistik. Pogoj za Neyman-Pearsonovo lemo je, da sta ničelna in alternativna domneva enostavni, torej natanko določata porazdeli- tev. 4.2 Moč testa Iz literature lahko povzamemo, da se športnikovo povprečje hemoglobina ob vsaj 14-dnevnem bivanju na višini nad 1500 m zviša za 2 g/l, medtem ko višinski treningi ne vplivajo na varianco njegovih vrednosti. Ob običajnih tre-ningih se posameznikove vrednosti porazdeljujejo normalno, X ∼ N (µ1, 52), kjer je µ1 športnikovo povprečje. Športnik pogosto opravlja višinske treninge, vendar v krajših intervalih. Za- nima ga, ali se njegovo povprečje hemoglobina v obdobju višinskih treningov kljub temu zviša. V sezoni opravi 12 meritev, 8 med obdobjem višinskih priprav in 4 sicer. Cilj naloge je ugotoviti, kolikšna bo moč njegovega testa, če bo pri sklepanju uporabil stopnjo značilnosti α = 0,05? a) Kaj je športnikova ničelna in kaj alternativna domneva? Zapišimo porazdelitev hemoglobina v obdobju višinskih priprav z N (µ2,52). Ničelna domneva je: H0: povprečje hemoglobina v obeh obdobjih je enako, µ1 = µ2. Alternativna domneva je, da je µ2 > µ1, zanima ga torej le enostranski test. b) Predlagajte testno statistiko. Izračunajte njeno porazdelitev pod ničelno domnevo. 78 Verjetnost in statistika 4.2. Moč testa Športnik bo izračunal razliko povprečij na vzorcih, ki je porazdeljena kot σ2 σ2 R = ¯ X2 − ¯ X1 ∼ N µ2 − µ1, + n1 n2 Testna statistika ¯ X Z = 2 − ¯ X1 σ2 + σ2 n1 n2 je torej pod ničelno domnevo porazdeljena standardno normalno. Ker ga zanima le enostranska alternativna domneva, bo ničelno domnevo zavrnil, kadar bo Z > zα, torej Z > 1,64. c) Izračunajte moč testa, torej verjetnost, da bo ničelno domnevo uspel za- vrniti, če se mu povprečje hemoglobina v obdobju višinskih priprav zares poveča za 2 g/l? Zanima nas PA(Z > 1,64). Označimo SE = σ 1 + 1 , pod alterna- n1 n2 tivno domnevo je pričakovana vrednost Z enaka   ¯ X µ E 2 − ¯ X1 2 − µ1 A(Z ) = E   = σ2 + σ2 SE n1 n2 torej velja Z ∼ N ( 2 ,1), in zato SE 2 PA(Z > 1,64) = P (U > 1,64 − ) SE kjer je U standardna normalna spremenljivka. V našem primeru   2 P U > 1,64 −  = P (U > 0,99) = 0,16 (4.1) 5 1 + 1 8 4 Moč testa je zelo majhna - pri tako majhnem številu meritev je le majhna verjetnost, da bo športnik zavrnil ničelno domnevo (četudi se mu pov- prečje dejansko zares zviša za 2 g/l). 79 Verjetnost in statistika 4.2. Moč testa d) Kako bi se moč testa spremenila, če bi imel na voljo enako število meritev v vsakem obdobju? Če bi imel po 6 meritev v vsakem obdobju, bi bila moč enaka   2 P U > 1,64 −  = P (U > 0,95) = 0,17 5 1 + 1 6 6 e) Kako je moč testa odvisna od variance posameznikovih meritev in kako od dejanske velikosti razlike v populaciji? Iz enačbe (4.1) je očitno, da večja razlika pomeni večjo moč - če je dejan- ska razlika med obdobjema večja, jo bomo lažje opazili na podatkih. Če bi bila varianca posameznikovih meritev manjša, bi imeli manjšo stan- dardno napako, in zato večjo moč testa. Povzetek • V nalogi smo si ogledali izračun moči testa pri primerjavi povprečij dveh skupin. Cilj opisane raziskave je dokazati, da obstaja razlika med skupinama, a smo z izračunom moči ugotovili, da je izvedba opisane raziskave precej nesmiselna, saj ima športnik le zelo majhno verjetnost, da bo razliko tudi dejansko dokazal. Za zanesljive trditve o razlikah bi potreboval precej večji vzorec. • Izračun moči testa je pomemben korak pri načrtovanju raziskave in ključ do izbire primerno velikega vzorca. Ponavadi pravimo, da razi- skave ni smiselno izvesti, če moč testa ni vsaj 0,8, torej če raziskova- lec nima vsaj verjetnosti 0,8, da bo ničelno domnevo uspel zavrniti. Za izračun moči moramo vedno nekaj vedeti o razpršenosti podatkov, hkrati pa moramo vedeti tudi, kako veliko razliko bi radi dokazali. To vprašanje ni vprašanje za statistika, temveč za raziskovalca - že pri načrtovanju raziskave mora vedeti, kako velike razlike bi rad dokazal oziroma kolikšna je najmanjša razlika, ki bo zanj tudi strokovno po- membna. • Če smo na podlagi vzorca s tveganjem, manjšim od α uspeli dokazati, da obstaja razlika med skupinama v populaciji, pravimo, da je razlika 80 Verjetnost in statistika 4.3. Enostavni domnevi statistično značilna. Vendar pa statistična značilnost še ne pomeni, da je razlika tudi dejansko strokovno pomembna oziroma da je smiselno razmišljati o razlogih zanjo. V populaciji bosta dve skupini le redko imeli natanko enako povprečje, zato bomo z zelo velikimi vzorci lahko skoraj vedno dobili statistično značilne rezultate. Kot primer si za- mislimo raziskavo, s katero želimo po višini primerjati prebivalce dveh mest. Povprečji obeh populacij zagotovo nista natanko enaki, lahko se na primer razlikujeta za 2 mm, kar pa je povsem nepomembno. Če bomo imeli zares velik vzorec, bomo na vzorcih torej lahko dobili sta- tistično značilno razliko, ki pa bo v praksi povsem nepomembna. • S statističnim testom nikoli ne bomo mogli pokazati, da v populaciji ni razlik med skupinama. A če navkljub veliki moči testa ne bomo uspeli zavrniti ničelne domneve, bo dobljeni interval zaupanja za razliko verjetno dovolj ozek, da bomo lahko morebitno razliko v populaciji omejili na vrednosti, ki niso strokovno pomembne. • Testna statistika, ki smo jo uporabili v nalogi, je zelo podobna stati- stiki, ki jo uporablja test t. Razlika je le, da v naši nalogi naredimo pre- cej predpostavk glede porazdelitve populacije (vemo, da je normalna, poznamo varianco), in zato lahko preprosto izračunamo tudi poraz- delitev testne statistike. Kadar teh predpostavk ne bomo naredili, si bomo do približne porazdelitve pomagali s centralnim limitnim izre- kom, upoštevati pa bomo morali tudi, da ne poznamo variance, temveč jo ocenjujemo s pomočjo vzorca. • V nalogi smo si ogledali še, kako na moč vpliva velikost skupin. Izkaže se, da je v smislu moči optimalen načrt raziskave tak, da sta skupini enako veliki. 4.3 Enostavni domnevi Računalnik skuša razlikovati med dvema viroma signalov. Prvi vir oddaja si- gnale, katerih jakost je normalno porazdeljena z N (0,1), drugi vir ima enako porazdelitev, a višjo povprečno jakost - N (2,1). Računalnik prejme 10 signalov in se mora odločiti, iz katerega vira so prišli. Posamezni signali iz istega vira so med seboj neodvisni. 81 Verjetnost in statistika 4.3. Enostavni domnevi a) Računalnik se odloča med domnevama H1: signal prihaja iz vira 1 in H2: signal prihaja iz vira 2 Zapišite testno statistiko, ob pomoči katere naj se odloča računalnik. (Iz- razite malo bolj splošno - naj imata porazdelitvi obeh virov varianco σ2, povprečna jakost drugega vira naj bo a, a > 0, velikost vzorca naj bo n.) Namig: pri tem, kaj je bolj verjetno, si pomagajte z gostotami. Kot testno statistiko uporabimo razmerje verjetij - bolj ko bo razmerje različno od 1, bolj bomo prepričani v eno izmed domnev: n n 1 f √ exp{− (xi−a)2 } 2(xi) 2σ2 = 2πσ f 1 1(xi) √ exp{− (xi)2 } i=1 i=1 2πσ 2σ2 n exp{− (xi−a)2 } = 2σ2 exp{− (xi)2 } i=1 2σ2 n (x = exp{− i − a)2 − (xi)2 } 2σ2 i=1 n −2ax = exp{− i + a2 } 2σ2 i=1 n 2ax = exp{ i − a2 } 2σ2 i=1 na a = exp{ ¯ x − } σ2 2 b) Denimo, da želimo, da računalnik reagira le, če je precej prepričan, da signal ne prihaja iz vira 1. Domnevo H1 proglasimo za ničelno domnevo, domnevo H2 pa za alternativno. Odločitveno pravilo postavimo tako, da bo verjetnost zmote, kadar je ničelna domneva resnična, največ α = 0,05. • Testna statistika je slučajna spremenljivka (označite jo z Y ). Kaj lahko rečemo o njeni porazdelitvi pod ničelno domnevo? Namig: da bo porazdelitev preprostejša, uporabite logaritem. Označimo Y = na ¯ X − a . Večja kot bo vrednost Y , bolj bomo σ2 2 prepričani v domnevo H2. Da bi vedeli, katere vrednosti so ‘velike’, moramo poznati porazdelitev slučajne spremenljivke Y . 82 Verjetnost in statistika 4.3. Enostavni domnevi Pod ničelno domnevo so vrednosti Xi standardno normalno poraz- deljene. Ker so a in σ2 konstante (znane vrednosti), je Y linearna kombinacija neodvisnih normalnih spremenljivk, in zato normalno porazdeljena. Poiščemo povprečje in standardni odklon spremen- ljivke Y : na na2 E ¯ 0(Y ) = E0 X − σ2 2σ2 na na2 na2 = E0( ¯ X) − = − σ2 2σ2 2σ2 na na2 var ¯ 0(Y ) = var0 X − = σ2 2σ2 n2a2 n2a2σ2 a2 = var ¯ 0X = = n σ4 σ4n σ2 √ a n sd0(Y ) = σ • Izrazite mejno vrednost, pri kateri naj računalnik reagira. Ničelno domnevo bomo zavrnili pri velikih vrednostih Y , zanima nas torej vrednost c, tako da bo P0(Y ≥ c) = 0,05, pri čemer P0 označuje verjetnost pri predpostavki, da ničelna domneva drži. Y je normalno porazdeljena slučajna spremenljivka, zato velja α = P0( Y −E(Y ) ≥ z sd(Y ) α) = α oz. P0(Y ≥ E(Y ) + zαsd(Y )). Za primer, ko je √ α = 0,05, je mejna vrednost c = − na2 +1,64 a n . Če je n = 10, a = 2 2σ2 σ √ in σ = 1, je mejna vrednost za c = − 10·4 + 1,64 · 2 10 = −9,63. 2 c) Kolikšna je verjetnost, da bo računalnik reagiral, če signal v resnici prihaja iz drugega vira? (Tej verjetnosti pravimo moč testa.) Namig: izračunajte porazdelitev testne statistike pod alternativno do- mnevo. Pod alternativno domnevo se Y prav tako porazdeljuje normalno (linearna kombinacija normalnih), povprečje je enako na na2 E ¯ A(Y ) = EA X − σ2 2σ2 na na2 na2 = a − = σ2 2σ2 σ2a 83 Verjetnost in statistika 4.3. Enostavni domnevi Ker je varianca pri obeh hipotezah enaka, je enaka tudi varianca vzorčnega povprečja, in zato je tudi varianca pod alternativno domnevo enaka a2 varA(Y ) = n σ2 Naj bo Z standardno normalno porazdeljena spremenljivka. Moč testa je √ na2 a n PA(Y > c) = PA(Y > − + zα ) 2σ2 σ √ Y − na2 − na2 + z a n α − na2 = P 2σ2 2σ2 σ 2σ2 A( √ > √ ) a n a n σ σ √ zα n − na = P (Z > σ √ ) n √na = P (Z > zα − ) σ Vidimo, da je moč testa odvisna od velikosti vzorca (večji vzorec, večja moč), variance (če podatki bolj variirajo, imamo manjšo moč) ter seveda od povprečja a pod alternativno domnevo. V našem primeru je povprečje pod alternativno domnevo kar za dva standardna odklona proč od pov- prečja pod ničelno domnevo, zato je moč zelo velika navkljub majhnemu vzorcu: √na √ P (Z > zα − ) = P (Z > 1,64 − 2 10) = P (Z > −4,68) σ Moč testa je skoraj enaka 1. Vidimo tudi, da je spodnja meja moči (a je večji od 0) enaka α, kar je verjetnost, da zavrnemo ničelno domnevo, kadar je a = 0. d) Testno statistiko Y transformirajte tako, da bo pod ničelno domnevo stan- darna normalna spremenljivka. 84 Verjetnost in statistika 4.3. Enostavni domnevi Standardizirajmo Y : Y + na2 Z = 2σ2 √ a n σ na ¯ X − na2 + na2 = σ2 2σ2 2σ2 √ a n σ ¯ X = √ σ/ n e) Povzemite: ali je mejna vrednost testne statistike odvisna od a? Intui- tivno razložite. Je vrednost a torej sploh pomembna? V prejšnji točki smo testno statistiko standardizirali, ker gre za bijektivno preslikavo (linearno transformacijo), je Z testna statistika, ki je ekviva-lentna Y (oz. prvotno definirani testni statistiki). S tem smo pokazali, da mejna vrednost testne statistike ni odvisna od vrednosti a. To je intui- tivno smiselno - mejna vrednost pod ničelno domnevo je postavljena glede na porazdelitev pod ničelno domnevo - poznavanje alternativne domneve za določitev mejne vrednosti ni potrebno (je pa pomembno, da vemo, da je enostranska). So pa vrednosti pod alternativno domnevo seveda ključne za moč testa. Predlogi za vaje v R • Generirajte podatke v dveh korakih. Najprej z enako verjetnostjo iz- berite vrednost a (0 ali 2), nato generirajte 10 vrednosti iz porazdeli- tve N (a,1). Izračunajte vrednost testne statistike iz prve točke in se odločite za eno izmed domnev glede na to, ali je vrednost testne stati- stike večja ali manjša od 1. Postopek velikokrat ponovite in izračunajte delež poskusov, v katerih se odločite za vsako od domnev, ter delež po- skusov, ko je ta odločitev pravilna. • Zamenjajte verjetnost v prvem koraku (naj bo npr. bolj verjetna izbira a = 2). Kako se spreminjajo deleži iz prejšnje točke? • Generirajte podatke, tako da je a ves čas enak 0, in preverite, da je vrednost α pri izračunani mejni vrednosti c zares enaka 0,05. • Generirajte podatke še tako, da je a = 2, in preverite moč testa. 85 Verjetnost in statistika 4.3. Enostavni domnevi Povzetek • V nalogi smo se želeli odločati med dvema domnevama o populaciji, v ta namen smo najprej definirali testno statistiko, ki povzame ustre- zno informacijo iz vzorca in katere vrednosti kažejo na eno ali drugo domnevo. Glede na to katera od domnev drži v populaciji, obstaja določena verjetnost, da testna statistika zavzame vrednost, ki smo jo opazili. Kriterij za odločanje med domnevama temelji na teh verjetno- stih. V tej nalogi eno izmed domnev določimo za našo osnovno, za drugo se odločimo le, če obstaja le majhna verjetnost, da vzorec izhaja iz populacije, v kateri drži osnovna domneva. • V nalogi smo izpeljali porazdelitev testne statistike na vzorcih iz popu- lacije, v kateri drži ničelna domneva. Na podlagi te porazdelitve smo določili kriterij za odločanje glede na željeno napako prve vrste. Na- sprotno nam porazdelitev testne statistike pod alternativno domnevo omogoča izračun moči našega testa. • V nalogi smo pokazali, da je moč testa odvisna od: – velikosti vzorca: večji kot je naš vzorec, bolj smo prepričani v svoje zaključke in manjša je verjetnost, da zagrešimo napako; – velikosti razlike v populaciji: večja bo razlika, lažje bomo razločili med skupinama tudi že pri manjših vzorcih; – variance v populaciji: če so razlike med enotami v populaciji majhne, tudi z majhnim vzorcem ne moremo pretirano zgrešiti prave ocene povprečja. Če pa je razpršenost v populaciji velika, bodo tudi vzorčne vrednosti zelo razpršene, naše ocene precej ne- natančne, in zato naši zaključki manj gotovi. • Opisani test imenujemo test razmerja verjetij. Seveda bi se lahko v danem primeru domislili tudi kake druge testne statistike, a uporabljena ima zelo pomembno lastnost. Neyman-Pearsonova lema (Rice, 2009, stran 332) namreč pove, da ima test razmerja verjetij pri odločanju med dvema enostavnima domnevama največjo moč med vsemi možnimi testi, ki bi jih lahko definirali pri enaki velikosti α. 86 Verjetnost in statistika 4.4. Enostavni domnevi, posplošitev 4.4 Enostavni domnevi, posplošitev Prejšnjo nalogo zapišimo v splošnem (Perman, 2013). Predpostavljamo, da so opazovane vrednosti neodvisne, enako porazdeljene slučajne spremenljivke X1, X2, ..., Xn. Predpostavite, da sta samo dve mož- nosti: ali je gostota spremenljivk enaka f (x) ali pa g(x), kjer sta f (x) in g(x) znani pozitivni gostoti. Formalno postavimo: H0: gostota je f (x) in H1: gostota je g(x) a) Predlagajte testno statistiko za preizkušanje zgornje domneve, če imate opazovane vrednosti x1, ..., xn. Kot testno statistiko uporabimo razmerje n g(X L = i) f (Xi) i=1 Velike vrednosti L so ‘dokazno gradivo’ proti ničelni domnevi. b) Kdaj bi zavrnili ničelno domnevo pri stopnji značilnosti α? Izrazite apro- ksimativno kritično mejo s količinama g(x) g(x) 2 a = log f (x)dx in b = log f (x)dx f (x) f (x) R R Zanima nas porazdelitev testne statistike pod ničelno domnevo, za testno statistiko vzemimo n g(X W = log i) f (Xi) i=1 Opazimo, da je testna statistika vsota slučajnih spremenljivk Yi = log g(Xi) , f (Xi) ki so neodvisne in enako porazdeljene, saj so take tudi spremenljivke Xi. Zato lahko uporabimo centralni limitni izrek - spremenljivka W je približno normalno porazdeljena (ne glede na porazdelitev Xi). Da bi lahko izračunali poljubno verjetnost, moramo poznati parametra te nor- malne porazdelitve. Pod ničelno domnevo lahko pričakovano vrednost Yi izračunamo kot g(x) E0(Yi) = log f (x)dx f (x) R 87 Verjetnost in statistika 4.5. Posplošeni test razmerja verjetij Pričakovana vrednost spremenljivke W je torej na. Podobno lahko izpe- ljemo, da je varianca spremenljivke W pod ničelno domnevo enaka var0(W ) = nvar0(Yi) = n(b − a2) Aproksimativno torej velja P (l > na + zα n(b − a2)) ≈ α. Povzetek • Naloga predstavlja posplošitev prejšnje - medtem ko sta v prejšnji na- logi obe domnevi predpostavljali normalno porazdelitev, smo tu vzeli povsem splošni porazdelitvi. V obeh primerih je razmerje gostot smi- selna testna statistika, vendar pa v splošnem njene porazdelitve seveda ne poznamo. • V nalogi smo približek za porazdelitev testne statistike poiskali s pomočjo centralnega limitnega izreka. Testno statistiko smo zapisali kot vsoto neodvisnih enako porazdeljenih slučajnih spremenljivk (funkcij slučajnih spremenljivk), zato konvergira proti normalni porazdelitvi. Izpeljali smo tudi parametra te normalne porazdelitve (povprečje in std. od- klon) in tako lahko določili ustrezne meje za testno statistiko. • Za nekatere porazdelitve f in g (npr. normalno, eksponentno ipd.) lahko izpeljemo natančno porazdelitev testne statistike. Kadar to ni mogoče, si pomagamo s približkom, ki velja asimptotsko. Vendar pa to ne zagotavlja, da bo test zadovoljivo deloval na majhnih vzorcih in pod ničelno domnevo dejansko zavračal z želeno napako α. Kadarkoli imamo torej le asimptotske približke testne statistike, moramo veljav- nost testa na majhnih vzorcih preveriti s simulacijami. 4.5 Posplošeni test razmerja verjetij Zanima nas, ali imajo zares vsi športniki enako variabilnost hemoglobina. Primerjati želimo meritve k športnikov, naj bodo vrednosti i-tega športnika (i = 1, . . . , k) porazdeljene normalno, torej Xij ∼ N (µi, σ2), kjer j = 1, . . . , n i i označujejo meritve pri posamezniku. Predpostavimo, da so vse meritve med seboj neodvisne. a) Zapišite ničelno in alternativno domnevo ter testno statistiko (posplošeni test razmerja verjetij). 88 Verjetnost in statistika 4.5. Posplošeni test razmerja verjetij Ničelna domneva: H0: σ2 = σ2 = . . . = σ2 1 2 k Alternativna domneva: H1 : σ2 niso vse enake i Testna statistika posplošenega testa razmerja verjetij je enaka: max L(x, θ) θ∈ω L(x, θ Λ = A∪ω0 = A∪0) max L(x, θ) L(x, θ θ∈ω 0) 0 kjer sta ω0 in ωA prostora parametrov, kot jih določa posamezna domneva, θ0 in θA pa oceni po metodi največjega verjetja glede na posamezno dom- nevo. Ker se števec podatkom lahko bolj ’prilagodi’ kot imenovalec, bo vrednost ulomka vedno večja od 1. Večja kot je, bolj je to dokaz v prid alternativne domneve, zanima nas, ali je vrednost Λ pod ničelno domnevo nenavadno velika. Da bi preverili ničelno domnevo, moramo torej narediti naslednje: zapisati funkcijo verjetja pod ničelno in alternativno domnevo, oceniti parametre po metodi največjega verjetja za obe domnevi, vstaviti cenilke za parametre v izraz za Λ in izračunati izraz na svojih podatkih. b) Najprej vzemimo, da imamo le enega športnika in n njegovih meritev. Kako bi ocenili njegova parametra µ in σ2 z metodo največjega verjetja? Funkcija verjetja je enaka n 1 (x L(x, µ, σ) = √ exp{− j − µ)2 } 2πσ 2σ2 j=1 del njenega logaritma, v katerem nastopata parametra, ki ju želimo oce- niti, pa je enak n 1 log L(x, µ, σ) = −n log σ − (xj − µ)2 2σ2 j=1 89 Verjetnost in statistika 4.5. Posplošeni test razmerja verjetij Poiščimo maksimum po µ: ∂ log L(x, ˆ µ, ˆ σ) = 0 ∂µ n 1 − (xj − µ)(−2) = 0 2σ2 j=1 n (xj − µ) = 0 j=1 n 1 µ = xj n j=1 Pa še za varianco: ∂ log L(x, ˆ µ, ˆ σ) = 0 ∂σ n n 1 −2 − − (xj − µ)2 = 0 σ 2 σ3 j=1 n −σ2n + (xj − µ)2 = 0 j=1 n 1 σ2 = (xj − µ)2 n j=1 c) Vrnimo se h k športnikom. Utemeljite, da so, kadar prostora možnih parametrov ne omejujemo z ničelno domnevo, ocene parametrov enake n 1 i ˆ µi = xij ni j=1 n 1 i ˆ σ2 = (x i ij − ˆ µi)2 ni j=1 Funkcija verjetja je enaka k ni 1 (x L(x, µ, σ) = √ exp − ij − µi)2 2πσ 2σ2 i=1 j=1 i i 90 Verjetnost in statistika 4.5. Posplošeni test razmerja verjetij Vsak člen vsote, ki jo dobimo po logaritmiranju gornje funkcije, je ses- tavljen le iz parametrov enega posameznika, ko odvajamo po tistem pa- rametru, ostanejo torej le členi, ki so vezani na tistega posameznika. Za ocenjevanje parametrov za nekega posameznika i torej potrebujemo iz- ključno njegove vrednosti, parametre posameznikov ocenimo povsem ne- odvisno drug od drugega. d) Kakšna je ocena povprečij pod ničelno domnevo? Pod ničelno domnevo je σi enak za vse i, zato ga v logaritmu funkcije verjetja lahko izpostavimo in ne vpliva na našo oceno posameznih pov- prečij. Ocena posameznih povprečij je zato enaka kot, če dovolimo različne variance. e) Kakšna je ocena variance pod ničelno domnevo? Del logaritma funkcije verjetja, ki nas zanima, je enak k k n 1 i log L(x, µ, σ0) = − ni log σ0 − (xij − µi)2. 2σ2 i=1 0 i=1 j=1 Odvod po σ izenačimo z 0 in dobimo k n 1 i ˆ σ2 = (x 0 ij − ˆ µi)2 k n i=1 j=1 i i=1 f) Kako bi ničelno domnevo preverili s testom razmerja verjetij? Zapišemo Wilksov Λ (zgoraj je funkcija verjetja na celotnem prostoru 91 Verjetnost in statistika 4.5. Posplošeni test razmerja verjetij parametrov, spodaj pod ničelno domnevo): k ni 1 √ exp{− (xij−µi)2 } 2πσ 2σ2 i=1 j=1 i i Λ = k ni 1 √ exp{− (xij−µ0i)2 } 2πσ0 2σ2 i=1 j=1 0 ni k n (x i k ij −µi)2 1 √ exp{− j=1 } 2πσ 2σ2 i=1 j=1 i i=1 i = k ni k ni 1 (x √ exp{− ij −µ0i)2 } 2πσ0 2σ2 i=1 j=1 i=1 j=1 0 Vstavimo ocene za variance v eksponent in tako v števcu kot v imenovalcu k dobimo exp{− 1 n 2 i}, ki se zato pokrajša. Logaritem Λ je enak i=1 k ni k ni log Λ = − log(σi) + log(σ0) i=1 j=1 i=1 j=1 k k = ni log(σ0) − ni log(σi) i=1 i=1 k = ni [log(σ0) − log(σi)] i=1 Dvakratna vrednost logaritma verjetij je porazdeljena kot χ2 , saj smo k−1 pod ničelno domnevo ocenili k − 1 parametrov manj. Predlogi za vaje v R • Generirajte primer s k športniki in ni meritvami pri vsakem športniku. Individualna povprečja športnikov naj bodo normalno porazdeljena okrog 148 s standardnim odklonom 7. Vzemite, da ničelna domneva drži in da so standardni odkloni okrog individualnih povprečij za vse športnike enaki (npr. 5). Na vsakem generiranem vzorcu ocenite po- trebne parametre in izračunajte vrednost testne statistike. Porazdelitev vzorčnih vrednosti testne statistike primerjajte s teoretično (asimptot- sko) porazdelitvijo. 92 Verjetnost in statistika 4.5. Posplošeni test razmerja verjetij • Pri delu si lahko pomagate s spodnjo kodo. Oglejte si, kako spremi- njanje parametrov vpliva na to, ali bo asimptotska porazdelitev upo- rabna pri vašem primeru - predvsem spreminjajte velikost vzorca, torej število športnikov in število meritev na športnika. Opazili boste, da je pri majhnih vzorcih asimptotska porazdelitev dokaj slab približek. Ali to v vašem primeru pomeni, da se vrednost α poveča ali zmanjša? • Spremenite še vrednosti varianc in izračunajte moč testa za nekaj pri- merov. > runs <- 100 #stevilo simulacij > rez<- rep(NA,runs) #pripravimo za rezultate > k <- 30 #stevilo sportnikov > ni <- 10 #stevilo meritev na sportnika > mi <- rnorm(k,148,7) #povprecna vrednost za vsakega sportnika > si <- rep(5,k) #variance za vsakega sportnika > for(jt in 1:runs){ #simulacija po korakih > x <- NULL #pripravimo za vpisovanje vrednosti > for(it in 1:k){ #za vsakega sportnika posebej > xij <- rnorm(ni,mi[it],si[it]) #generiramo vrednosti > x <- c(x,xij) #zdruzimi s prejsnjimi > } > data <- data.frame(id=rep(1:k,each=ni),x=x) #podatki #sedaj pa ocenimo parametre s pomocjo podatkov > mihat <- sihat <- rep(NA,k) > for(it in 1:k){ #za vsakega sportnika > mihat[it] <- 1/ni*sum(data$x[data$id==it]) #ocena povprecja #ocena standardnega odklona: > sihat[it] <- sqrt(1/ni * sum((data$x[data$id==it] - mihat[it])^2)) > } > mihat <- rep(mihat,each=ni) #ocena standardnega odklona pod nicelno domnevo > s0hat <- sqrt(1/(ni*k)* sum((data$x - mihat)^2)) #log Lambda: > logL <- ni*sum(log(s0hat) - log(sihat)) > rez[jt] <- 2*logL #vrednost testne statistike na tem vzorcu > } 93 Verjetnost in statistika 4.5. Posplošeni test razmerja verjetij > plot(ecdf(rez)) #empiricna porazdelitvena funkcija > x <- seq(min(rez),max(rez),length=100) > lines(x,pchisq(x,k-1),col=4) #teoreticna porazdelitvena funkcija Povzetek • Posplošeni test razmerja verjetij predstavlja generičen pristop, ki nam poda smiselno testno statistiko in njeno porazdelitev. Testna statistika je enaka kvocientu maksimumov dveh verjetij - na celotnem prostoru parametrov in na podprostoru, ki ga določa ničelna domneva. V našem primeru je celoten prostor parametrov dovoljeval k vrednosti za vari- anco, medtem ko smo pod ničelno domnevo zahtevali, da so vse te vrednosti enake. Nasprotno smo v primeru iz naloge 2.2 v števcu vzeli le vrednost pod alternativno domnevo in ne vseh možnih vrednosti pa- rametra (torej pod ničelno in pod alternativno domnevo). • Porazdelitev testne statistike pri posplošenem testu razmerja verjetij je seveda le približna in na majhnih vzorcih je lahko približek precej slab. To pa pomeni, da bomo ob uporabi te porazdelitve na majhnih vzorcih zavračali z verjetnostjo, ki ne bo povsem enaka α. To smo raziskali s primerom v R in pokazali, da na majhnih vzorcih (pri majhnem številu meritev na športnika) ničelno domnevo zavračamo prepogosto. • Mimogrede opazimo tudi, da ocena σ po metodi največjega verjetja ni nepristranska - metoda največjega verjetja zagotavlja le doslednost (cenilka konvergira k pravi vrednosti). 94 Poglavje 5 Linearna regresija To poglavje predstavlja le kratek uvod v linearno regresijo (Rice, 2009, po- glavje 14). V prvi nalogi je preizkušanje domnev v linearni regresiji predstavljeno le kot poseben primer, pri katerem lahko uporabimo različne testne statistike z različnimi lastnostmi. Ogledali si bomo Waldov test ter izpe- ljali testno stastistiko za posplošeni test razmerja verjetij. Druga naloga je posvečena matričnemu računanju ter izpeljavi cenilk in njihovih lastnosti v linearni regresiji. Tovrstni pristop je neprimerno bolj pregleden od klasičnega, in zato nujen za vsakršne zahtevnejše izpeljave. Za konec sledi še naloga, ki se dotakne predpostavk linearne regresije in nakaže težave, ki se pojavijo, če so posamezne predpostavke kršene. 5.1 Linearna regresija Zanima nas povezanost števila ur učenja na teden z rezultatom na izpitu iz statistike. Vzemimo, da vemo, da se rezultat na izpitu v populaciji porazde- ljuje pogojno normalno: Y |X ∼ N (β0 + β1X, σ2). a) Iz spodnjega izpisa preberite ocene populacijskih parametrov. Interpreti- rajte rezultate, katere ničelne domneve so testirane in kako? Call: lm(formula = y ~ x) Residuals: Min 1Q Median 3Q Max -20.683 -4.746 2.844 4.512 14.693 95 Verjetnost in statistika 5.1. Linearna regresija Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 19.2049 7.5172 2.555 0.033921 * x 3.6850 0.6217 5.927 0.000351 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 11.44 on 8 degrees of freedom Multiple R-squared: 0.8145, Adjusted R-squared: 0.7913 F-statistic: 35.13 on 1 and 8 DF, p-value: 0.0003508 Ocene parametrov so β0 = 19,2, β1 = 3,7, σ = 11,4. Testirani sta dve ničelni domnevi: H0int: β0 = 0 in H0: β1 = 0. Pri linearni regresiji nas ponavadi zanima le druga - saj ta govori o povezanosti med spremenljiv- kama v populaciji. Pri iskanju porazdelitve cenilke β1 bi se lahko oprli na teorijo metode največjega verjetja, vendar pa v tem primeru aproksi- macija ni potrebna. Cenilka je namreč linearna kombinacija vrednosti Y (glej nalogo 3.2), zato je normalno porazdeljena. Njena varianca (stan- dardna napaka) je ocenjena iz podatkov, zato je standardizirana vrednost cenilke porazdeljena kot t. Testna statistika β 3,7 T = 1 = = 5,9 SE 0,6 β1 je torej porazdeljena kot t z 8 stopinjami prostosti (pri ocenjevanju SE porabimo dve stopinji prostosti). Ta test se imenuje Waldov test. b) Kako bi ničelno domnevo H0: β1 = 0 preverili s posplošenim testom raz- merja verjetij? Namig: kjer je le mogoče, uporabite rezultate iz prejšnje naloge. Začnimo z ocenami pod ničelno domnevo. Pod ničelno domnevo je pov- prečje za vse posameznike enako, neposredno torej lahko uporabimo rezul- tate iz naloge 4.5, le da namesto µ pišemo β0, zato je maksimum funkcije 96 Verjetnost in statistika 5.1. Linearna regresija verjetja pod ničelno domnevo enak 1 n (yi − β0)2 L i=1 0(y, x, β0, σ) = √ exp − ( 2πσ)n 2σ2 1 n n (yi − β0)2 = √ exp − i=1 ( 2πσ)n 2 n (y i=1 i − β0)2 1 n = √ exp − ( 2πσ)n 2 Za števec bomo uporabili rezultat, da je ocena σ enaka n (yi − β0 − xiβ1)2 σ2 = i=1 n Funkcija verjetja je enaka: 1 n (yi − β0 − β1xi)2 L(y, x, β i=1 0, β1, σ) = √ exp{− } ( 2πσ)n 2σ2 in tako je maksimum funkcije verjetja enak 1 n (yi − β0 − β1xi)2 L i=1 A∪0(y, x, β0, β1, σ) = √ exp − ( 2πσ)n 2σ2 1 n n (yi − β0 − β1xi)2 = √ exp − i=1 ( 2πσ)n 2 n (y i=1 i − β0 − β1xi)2 1 n = √ exp − ( 2πσ)n 2 Wilksov Λ je enak 1 L √ exp − n ( 2πσ 2 Λ = A∪0 = A)n L 1 0 √ exp − n ( 2πσ0)n 2 σn = 0 σn A  n n (yi − β00)2 =  i=1   n   (y  i − β0A − xiβ1A)2 i=1 97 Verjetnost in statistika 5.1. Linearna regresija Vrednost maksimuma pri neomejenem prostoru parametrov izračunamo tako, da vstavimo ocenjene β0 in β1, za izračun vrednosti pod ničelno domnevo moramo oceniti še β0 v ničelnem modelu. Vrednost 2 log Λ se porazdeljuje kot χ2. 1 Predlogi za vaje v R • Naj bo X enakomerno porazdeljena spremenljivka (med 0 in 20, zao- krožena navzdol), β0 = 15, β1 = 4, σ = 10. Generirajte vzorec velikosti 10, narišite podatke in vrišite populacijsko ter ocenjeno vrednost pre- mice. > set.seed(1) > n <- 10 #velikost vzorca > beta0 <- 15 > beta1 <- 4 > sigma <- 10 > x <- floor(runif(n)*20) #navzdol zaokrozene vrednosti x > x <- sort(x) #uredimo podatke po velikosti x > y <- rnorm(n,mean=beta0+beta1*x,sd=sigma) #generiramo y-one > plot(x,y) #narisemo tocke > popul <- beta0 + beta1*x #populacijska vrednost premice > lines(x,popul,col="grey",lwd=2) #dodamo popul. vrednost premice > fit <- lm(y~x) #ocenimo premico na podatkih > summary(fit) #ogledamo si ocene koeficientov > beta0h <- fit$coef[1] #ocenjena beta0 > beta1h <- fit$coef[2] #ocenjena beta1 > napoved <- beta0h + beta1h*x > lines(x,napoved,lwd=2) #vrisemo ocenjeno premico na sliko • Izračunajte posplošeni test razmerja verjetij v R > fit0 <- lm(y~1) #pod nic. domnevo - le konstanta > res0 <- y - fit0$coef #ostanki pod nicelno domnevo > resA <- y - beta0h - beta1h*x #ostanki pod alternativno domnevo #zanima nas razlika log verjetij - konstanto lahko izpustimo: > logl0 <- -.5*n*log(sum(res0^2)) #loglik pod nicelno > loglA <- -.5*n*log(sum(resA^2)) #loglik pod alternativno > Lambda <- 2*(loglA-logl0) #Wilksov lambda > 1-pchisq(Lambda,1) #likelihood ratio test [1] 4.048e-05 98 Verjetnost in statistika 5.2. Matrično računanje Povzetek • V nalogi smo si ogledali osnovni primer linearne regresije z eno neod- visno spremenljivko. Najprej smo si ogledali izpis programa R in ga interpretirali. Omenili smo dve možnosti preizkušanja domneve o re- gresijskem koeficientu: Waldov test in posplošeni test razmerja verjetij. • Pri obeh omenjenih testih je znana le asimptotska porazdelitev testne statistike, ničesar pa nismo povedali o kvaliteti aproksimacije na majh- nih vzorcih. Izkaže se, da ima posplošeni test razmerja verjetij na majhnih vzorcih pogosto boljše lastnosti od Waldovega testa, ki je av- tomatično poročan v vseh statističnih programih. 5.2 Matrično računanje Vrednosti neodvisnih spremenljivk združimo v matriko X (design matrix), vrednosti odvisne spremenljivke ter koeficientov predstavljajo vektorja Y in β:       1 x11 . . . x1p Y1 β0  1 x21 . . . x2p   Y2   β1  X =  . . .  ; Y =  .  ; β =  .   .. .. . . . ..   ..   ..        1 xn1 . . . xnp Yn βp Matrika X je dimenzije n × (p + 1), kjer je p število spremenljivk. Če naš model ne bi vseboval konstante, bi prvi stolpec X izpustili. n a) Zapišite vsoto vrednosti Y 2 v matrični obliki. i i=1   Y1  Y2  Y T Y = Y   1 Y2 . . . Yn . = Y 2 + Y 2 + . . . + Y 2  . 1 2 n .    Yn b) Kaj dobimo, če matrično pomnožimo Xβ (da bo manj pisanja, vzemite p = 1)? 99 Verjetnost in statistika 5.2. Matrično računanje       1 x1 β0 + β1x1 E(Y1)  1 x2  β  β0 + β1x2   E(Y2)  Xβ = 0  . .  =  .  =  .  = E(Y )  .. ..  β  .   .  1 . .       1 xn β0 + β1xn E(Yn) c) V matrični obliki oceno koeficientov po metodi najmanjših kvadratov (= po metodi največjega verjetja) zapišemo kot β = (XT X)−1XT Y . Pokažite, da za p = 1 dobite oceni: n n n n xiyi − xi yi β i=1 i=1 i=1 1 = n n n x2 − ( x i i)2 i=1 i=1 n n n n x2 y x x i i − i iyi β i=1 i=1 i=1 i=1 0 = n n n x2 − ( x i i)2 i=1 i=1 Izračunajmo najprej XT X:   1 x1  n  n x 1 1 . . . 1 i  1 x2  XT X =  . .  =  i=1  x n n  . .    1 x2 . . . xn . .   xi x2i 1 xn i=1 i=1 Inverz te 2 × 2 matrike je enak:  n n  x2 − x 1 i i (XT X)−1 =  i=1 i=1  n n n   n x2 − ( x − x i i)2 i n i=1 i=1 i=1 Izračunajmo še XT Y :   Y1  n  Y 1 1 . . . 1 i  Y2  XT Y =  .  =  i=1  x n  .    1 x2 . . . xn .   xiYi Yn i=1 100 Verjetnost in statistika 5.2. Matrično računanje Velja torej  n n   n  x2 − x Y 1 i i i (XT X)−1XT Y =  i=1 i=1   i=1  n n n n     n x2 − ( x − x x i i)2 i n iYi i=1 i=1 i=1 i=1  n n n n  x2 Y x x 1 i i − i iYi =  i=1 i=1 i=1 i=1  n n n n n   n x2 − ( x − x Y x i i)2 i i + n iYi i=1 i=1 i=1 i=1 i=1 Zgornja vrstica pri tem predstavlja oceno β0, spodnja pa β1. d) Izpeljite oceno po metodi najmanjših kvadratov še v matrični obliki. Pri tem boste potrebovali naslednje formule za matrično računanje: (A + B)T = AT + BT ; (AT )T = A; (AB)T = BT AT ; ∂βT A ∂βT AT Aβ = A; = 2AT Aβ ∂β ∂β Namig: kaj minimiziramo? Kako zapišemo vsoto kvadriranih ostankov v matrični obliki? Če je v modelu ena spremenljivka, iščemo minimum funkcije n (Yi − β0 − β1xi)2 i=1 V matrični obliki iščemo minimum funkcije (Y − Xβ)T (Y − Xβ) = (Y T − βT XT )(Y − Xβ) = Y T Y − βT XT Y − Y T Xβ + βT XT Xβ = Y T Y − 2βT XT Y + βT XT Xβ Pri tem smo v zadnji vrstici uporabili, da je (βT XT Y )T = βT XT Y , saj je matrika dimenzije 1 × 1. Sedaj odvajamo po β ∂ (Y TY − 2βTXTY + βTXTXβ) = −2XTY + 2XTXβ ∂β 101 Verjetnost in statistika 5.2. Matrično računanje in izenačimo z 0 (ter predpostavimo, da XT X ni singularna): −2XT Y + 2XT Xβ = 0 XT Xβ = XT Y β = (XT X)−1XT Y e) Pokažite, da je ocena koeficientov nepristranska (vzemite, da so vrednosti x-ov dane in ne slučajne). E(β) = (XT X)−1XT E(Y ) = (XT X)−1XT Xβ = β f) Izpeljite formulo za standardno napako ocenjenih koeficientov v matrični obliki. Intuitivno razložite, od česa je odvisna standardna napaka koefi- cienta β1 (za p = 1). Uporabite formulo: var(cY ) = c varY cT . var(β) = var[(XT X)−1XT Y ] = (XT X)−1XT varY [(XT X)−1XT ]T = σ2(XT X)−1XT X(XT X)−1 = σ2(XT X)−1 Izpeljali smo variančno-kovariančno matriko, variance so na diagonali. Standardna napaka SEβ je torej enaka 1 σ2 SEβ = 1 (XT X)−1 22 nσ2 = n n n x2 − ( x i i)2 i=1 i=1 σ2 = n (xi − ¯ x)2 i=1 σ2 σ = = √ nσ2 σ n x x Standardna napaka koeficienta je tako kot vedno odvisna od velikosti vzorca n ter razpršenosti podatkov. Vrednost σ je standardni odklon 102 Verjetnost in statistika 5.2. Matrično računanje ostankov okrog premice - večja kot je, bolj se lahko zmotimo pri oceni premice. Vendar tu ni pomembna le absolutna velikost variance ostan- kov, zanima nas variabilnost ostankov glede na variabilnost neodvisne spremenljivke. Če je razpon x-ov majhen, je naša ocena pri isti varia- bilnosti ostankov manj natančna. Razložimo to na našem primeru - če bi v vzorec zajeli le posameznike, ki so se učili 3-5 ur, bi bila poveza- nost med spremenljivkama (npr. merjena s korelacijskim koeficientom) pri istih regresijskih koeficientih dosti manjša, in zato bi bila možna večja odstopanja pri ocenjevanju. g) Kako bi izračunali interval zaupanja za napovedano premico v našem pri- meru (p = 1)? Kako bo tak interval izgledal na sliki? Izračunamo standardno napako za vsako točko posebej. var(yi) = var(β0 + β0xi) = var(β0) + x2var(β i 1) + 2xicov(β0, β1) Matrični izračun (za vse točke naenkrat): var(Y ) = var(Xβ) = Xvar(β)XT = σ2X(XT X)−1XT ● 100 ● 80 ● ● ● 60 ● Y ● ● 40 ● 20 ● 5 10 15 X Slika 5.1: Točke na vzorcu in ocenjena premica z intervalom zaupanja. h) Recimo, da nas zanima, kako sta število ur učenja in spol (0 = ženski, 1 = moški) povezana z rezultatom na izpitu iz statistike. Predpostavimo model, ki vključuje interakcijo. Kako bi preverili, ali je število ur učenja pri moških povezano z rezultatom na izpitu? 103 Verjetnost in statistika 5.2. Matrično računanje Model, ki ga prilagodimo podatkom, zapišemo kot: Y = β0 + β1spol + β2ure + β3(ure ∗ spol) Ničelna domneva, ki jo želimo preveriti, je torej H0 : β2 + β3 = 0. Zapišemo v matrični obliki. Naj bo vektor aT = [0,0,1,1], zanima nas H0: aT β = 0. Varianca aT β je enaka var(aT β) = aT varβa, za preizkušanje ničelne domneve uporabimo Waldov test. Predlogi za vaje v R • V R v matrični obliki zapišite podatke, izračunajte oceno po metodi najmanjših kvadratov ter jo primerjajte z izpisom funkcije lm. • Ocenite tudi standardno napako ter jo primerjajte z izpisom • Narišite sliko napovedane premice ter ji dodajte interval zaupanja. > X <- cbind(1,x) > sigma <- summary(fit)$sigma > inv <- solve(t(X)%*%X) > mat <- X%*%inv%*%t(X) > se <- sigma*sqrt(diag(mat)) > betah <- c(beta0h,beta1h) > plot(x,y) > lines(x,X%*%betah) > t8 <- qt(.975,8) > lines(x,X%*%betah - t8*se,col=2) > lines(x,X%*%betah + t8*se,col=2) Povzetek • Predstavitev podatkov v linearni regresiji z vektorji in matrikami se izkaže za zelo učinkovit ter pregleden način izpeljevanja cenilk ter do- kazovanja njihovih lastnosti. Poleg tega so izpeljave splošnejše, saj niso odvisne od števila neodvisnih spremenljivk. V tej nalogi smo si zato ogledali osnove matričnega računanja ter jih uporabili za izpeljavo cenilk in lastnosti v linearni regresiji. • Ogledali smo si tudi, kako lahko Waldov test zapišemo splošneje in z njim preverimo poljubno domnevo o koeficientih v modelu, ki jo lahko zapišemo kot linearno kombinacijo parametrov. • Program R podatke vedno obravnava kot vektorje oziroma matrike, zato lahko cenilke v matrični obliki tudi neposredno uporabimo v kodi. 104 Verjetnost in statistika 5.3. Predpostavke linearne regresije 5.3 Predpostavke linearne regresije Z osnovnim modelom linearne regresije naredimo štiri predpostavke: • ostanki so okrog premice porazdeljeni normalno • varianca ostankov ni odvisna od vrednosti neodvisne spremenljivke (ho- moskedastičnost) • ostanki so med seboj neodvisni • povezanost med X in Y je linearna Kaj se zgodi z ocenami koeficientov, njihovo pričakovano vrednostjo, standar- dno napako in intervali zaupanja, če je katera izmed prvih treh predpostavk kršena? a) Kaj se spremeni v izpeljavah, če ostanki okrog premice niso porazdeljeni normalno? V tem primeru ocena koeficientov po metodi največjega verjetja ni enaka oceni po metodi najmanjših kvadratov. Ocena po metodi najmanjših kvadratov bo identična kot do sedaj, enaka bo tudi ocena standardne napake. Prav tako bo cenilka po metodi najmanjših kvadratov nepri- stranska ocena populacijskih vrednosti. Da pa bi izračunali kak interval zaupanja, moramo izpeljati porazdelitev cenilke. Vemo, da je cenilka po metodi najmanjših kvadratov linearna kombinacija vrednosti Y , zato smo pri normalni porazdelitvi (pri danih x-ih) ostajali znotraj normalne po- razdelitve. V splošnem to seveda ni res. Druga možnost je, da izpeljemo raje oceno po metodi največjega verjetja - v tem primeru se cenilke spremenijo, vemo pa, da je njihova asimptotska porazdelitev normalna. b) Recimo, da je varianca ostankov odvisna od x. Če varianca ostankov ni enaka za vsak x, moramo varianco pisati v ma- 105 Verjetnost in statistika 5.3. Predpostavke linearne regresije trični obliki, npr.:   w1 0 . . . 0 0  0 w2 . . . 0 0  Σ = σ  . . . . .   .. .. .. .. ..    0 0 . . . 0 wn Varianca sicer ne vpliva na oceno koeficientov po metodi najmanjših kva- dratov, vendar pa se spremeni ocena po metodi največjega verjetja, saj moramo maksimizirati funkcijo −2(Y − Xβ)T Σ−1(Y − Xβ): (Y − Xβ)T Σ−1(Y − Xβ) = = (Y T − βT XT )Σ−1(Y − Xβ) = Y T Σ−1Y − βT XT Σ−1Y − Y T Σ−1Xβ + βT XT Σ−1Xβ = Y T Σ−1Y − 2βT XT Σ−1Y + βT XT Σ−1Xβ in zato −2XT Σ−1Y + 2XT Σ−1Xβ = 0 XT Σ−1Xβ = XT Σ−1Y β = (XT Σ−1X)−1XT Σ−1Y Ustrezno se spremeni tudi varianca ocene: var(β) = var[(XT Σ−1X)−1XT Σ−1Y ] = (XT Σ−1X)−1XT varY [(XT Σ−1X)−1XT ]T = σ2(XT Σ−1X)−1XT Σ−1X(XT Σ−1X)−1 = σ2(XT Σ−1X)−1 Če so vrednosti wi znane, se v ocenjevanje koeficientov in standardne napake le vrine diagonalna matrika. Statistično sklepanje je enako kot do sedaj. c) Kaj pa če ostanki med seboj niso neodvisni? Potem variančna matrika Σ ni več diagonalna (je npr. bločno diagonalna). Rezultati bodo podobni tistim v prejšnji točki, bo pa seveda ocenjevanje odvisno od tega, kaj vemo o elementih Σ. 106 Verjetnost in statistika 5.3. Predpostavke linearne regresije Povzetek • Naloga predstavlja le uvod v zelo široko področje posplošitev line- arne regresije. Osnovni primer linearne regresije namreč zahteva precej stroge predpostavke, ki so v praksi seveda pogosto kršene. V nalogi smo si le površno ogledali, kaj se v teoretičnih izpeljavah spremeni, kadar je kršena posamezna predpostavka, medtem ko se z možnimi rešitvami oziroma njihovimi lastnostmi nismo ukvarjali. 107 Literatura Blagus, R. (2011). Razvrščanje visoko-razsežnih neuravnoteženih podatkov : doktorsko delo. Medicinska fakulteta, Univerza v Ljubljani. Casella, G. and Berger, R. L. (1990). Statistical inference. Wadsworth, Belmont, CA. Perman, M. (2013). Metodologija statističnega raziskovanja, izpitne naloge. URL: http://valjhun.fmf.uni-lj.si/ mihael/ul/vs/izpiti.html. R Development Core Team (2013). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. Rice, J. A. (2009). Mathematical Statistics and Data Analysis. Duxbury advanced series, Cengage learning. Shao, J. (2003). Mathematical Statistics, second edition. Springer-Verlag, New York. 108 Stvarno kazalo Asimptotska aproksimacija, 35 Metoda najmanjših kvadratov, 70 binomske porazdelitve, 37–39 Metoda največjega verjetja, 64–73 Moč testa, 76, 78–81, 83 Cenilka, 41 dosledna, 64 Napak druge vrste, 74 nepristranska, 41–44, 49, 52–53 Napaka prve vrste, 74 Centralni limitni izrek, 35–40, 45, 87 Neyman-Pearsonova lema, 86 Delež Območje zavrnitve, 76 ocena, 65–69 Domneva Pokritje, 68 alternativna, 76, 82 Porazdelitev enostavna, 76 χ2, 4, 92 ničelna, 74, 75, 82 normalna, 2, 7 Bernoullijeva, 8 Fisherjeva informacija, 64, 67 binomska, 14, 37–39 matrika, 71 eksponentna, 10 enakomerna, 8 Generiranje slučajnih spremenljivk, 8– gama, 4, 16 11 standardna normalna, 3 vzorčnega povprečja, 22–24 Interakcija, 103 Pravila matričnega računanja, 101 Interval zaupanja, 44 Pričakovana vrednost, 24–30 Kovarianca in korelacija, 27–29 ocena, 42 ocenjevanje, 51–55 Problem večkratnega preizkušanja, 4, 7 Linearna regresija, 95 matrični pristop, 99–105 Razmerje verjetij, 78, 82 ocenjevanje koeficientov, 69–73 posplošeni test, 88–94 predpostavke, 105–107 Standardna napaka, 23, 42–43, 46 testi, 95–98 Stopnja značilnosti, 74, 76 109 Verjetnost in statistika Stvarno kazalo Testna statistika, 82 Transformacija integrala verjetnosti, 11 Varianca, 24–30 ocena, 43–44, 47 Vsota slučajnih spremenljivk diskretnih, 12–14 odvisnih, 21 zveznih, 14–18 Vzorčenje končna populacija, 45–48 neskončna populacija, 42–45 slučajno, 41 stratificirano, 48–50 večstopenjsko, 58–63 Waldov test, 96 Zbir, 64 110 Document Outline Kolofon_Maja Pohar Perme verstatfin2