MC meetodid, praktikum 3

Tööleht on valminud Hariduse Infotehnoloogia Sihtasutuse IT Akadeemia programmi toel

Teema 1: empiirilise jaotusfunktsiooni ja vaadeldavale jaotusele vastava jaotusfunktsiooni erinevuste mõõtmisel põhinevad testid

Empiiriline jaotusfunktsioon annab iga piirkonna \(A_x=(-\infty,x)\) korral sinna piirkonda sattunud vaatluste arvu jagatise kõikide vaatluste arvuga, \[ F_n(x) := \frac{\#\{i \mid x_i \leq x \}}{n}\] Juhusliku suuruse \(X\) jaotusfunktsioon aga annab iga \(x\) korral sellisesse piirkonda sattumise tõenäosus: \[F(x)=P(X\in A_x)\] Kolmogorov-Smirnovi test mõõdab maksimaalset erinevust, st leiab valimi põhjal \[ K_n := \sup_x |F_n(x) - F(x)|\] ning annab ka tõenäosuse, et valitud jaotusele vastava juhusliku suuruse puhul \(n\) katse tegemisel saame vähemalt sama suure väärtuse.

Cramer - von Mises test ehk \(\omega^2\)-test mõõdab empiirilise jaotusfunktsiooni ja jaotusele vastava jaotusfunktsiooni (mille tihedusfunktsiooniks olgu \(f\)) vahelist erinevust, arvutades allpool toodud integraali erinevuse ruudust: \[\omega_n^2 := n \int_{-\infty}^\infty (F(x) - F_n(x))^2 f(x) dx\] Jällegi on teada, kuidas see erinevus käitub juhul, kui vaadeldakse juhuslikule suurusele \(X\) vastavate sõltumatute katsete tulemusi.

Harjutus 1

Kontrolli Kolmogorov-Smirnovi testiga, kas käskudega set.seed(20230919);x <- -5*log(1-runif(200)) tekitatud pseudojuhuslikud arvud võivad vastata jaotusele \(Exp(0.2)\). Kontrolli samuti, kas need võivad vastata ka jaotusele \(Exp(0.19)\)

Vihjed: Kolmogorov-Smirnovi testi rakendamiseks on käsk ks.test, millele tuleb anda valim, vastava jaotuse jaotusfunktsiooni nimi R-is ning lisaks jaotusfunktsioonile täiendavad parameetrid (näiteks normaaljaotuse korral mean=1, sd=2).

Harjutus 2

Korda eelmise harjutuse tulemusi Cramer-von Mises testiga.

Vihjed: selle testi kasutamiseks peab R-is olema installeeritud lisapakett goftest, mis tuleb pärast installeerimist laadida. Testi sooritamise käsuks on cvm.test, mille kasutamin on sarnane Kolmogorov-Smirnovi testiga.

Teema 2: Sõltumatuse testimine

Eelnevalt vaadeldud testid kontrollisid ainult seda, kas juhusliku suuruse väärtused paiknevad võimalike tulemuste hulgas nii, nagu teoreetilise jaotuse järgi oodata võiks, kuid ei arvestanud seda, kuidas need väärtuse katseseerias esinesid. Seega võivad katsetulemuste vahelised sõltuvused avastamata jääda.

Siit tuleb idee vaadelda mitte ainult katsetulemusi ükshaaval, vaid mingite komplektidena ja vaadelda komplektidega seotud sündmuste vastavust sõltumatute katsete korral saadavatele tulemustele.

Esimene idee - sõltumatute paaride test. Selleks

Fikseerime \(m\) vastastikku välistavat sündmust \(A_1,\ldots, A_m\), mille tõenäosused vaadeldava jaotuse korral on \(p_1,\ldots,p_m\) ja mis katavad ära kõik võimalused (täissüsteem)
Igale katsetulemusele seame vastavusse sündmuse numbri, mis sellel katsel toimus
Saame numbritest \(1,..,m\) koosneva jada \(k_1,k_2,\ldots,k_n\)
Vaatleme paare \((k_1,k_2),(k_3,k_4),\ldots\)
Kui tegemist sõltumatute katsetega, siis paari \((i,j)\) saamise tõenäosus on \(p_i\cdot p_j\)
Kasutame \(\chi^2\) testi kontrollmiseks, kas esinemiste arvud \(n_{i,j}\) on piisavalt lähedased oodatavale arvule \(\frac{n}{2}p_i p_j\)
NB! testi rakendamiseks sagedustabeli koostamisel peame olema hoolikad, et kõik paarid tabelis esineksid (ka siis, kui katsetulemuste põhjal sellist paari ei tekkinud)! Samuti peaks veenduma, et iga sündmuste paari korral on toimumiste oodatavad arvud \(\frac{n}{2}p_i p_j\) piisavalt suured (vähemalt 5), vastasel korral ei ole õigustatud \(\chi^2\)-testi rakendamisel \(p\)-väärtuseid leida \(\chi^2\)-jaotuse abil vaid need tuleks leida näiteks simulatsioonide teel (suvand simulate.p.value=TRUE käsus chisq.test)

Harjutus 3

Rakenda sõltumatute paaride testi harjutuses nr 1 kasutatud jada järjestikuste vaatluste sõltumatuse kontrolliks sündmuste \(\{X\leq 2\},\{2<X\leq 4\},\{X>4\}\) korral (kasutades kontrollitava jaotusena \(Exp(0.2)\)).

Lahenduse skeem: Tekita toodud sündmustele vastav faktortunnus (olgu selle nimeks sildistatud_x) ja sündmuste tõenäosuste vektor (olgu selle nimeks p). Tekita faktorist sildistatud_x kaks faktorit esimene ja teine, kus esimene vastab paaride esimestele elementidele ja teine paaride teistele elementidele. Selleks tuleb faktorisse esimene korjata sildistatud_x esimene, kolmas jne element; faktorisse teine aga teine, neljas jne element. Nii vektoritest kui faktoritest kindlate elementide väljanoppimiseks võib indeksisse kirjutada vajaminevate koordinaatidega vektor ning praegusel juhul sobib selle tekitamiseks hästi käsk seq.

Seejärel saab lugeda moodustada tabeli, mille i-ndas reas ja j-ndas tulbas on paaride arv, kus esimene element vastas sündmusele \(A_i\) ja teine element sündmusele \(A_j\). Selleks sobib jälle table käsk, millele tuleb nüüd anda kaks argumenti (faktorid esimene ja teine).

Nüüd tee valmis tõenäosuste maatriks, mille element indeksiga \(i,j\) on \(p_i\cdot p_j\). Sellise tulemuse annab \(p\) komponentidest moodustatud üheveerulise maatriksi korrutamine \(p\) komponentidest moodustatud üherealise maatriksiga. Sarnast tegevust maatriksi tekitamiseks kahest vektorist läheb mitmel pool vaja ja seetõttu on vastaval operatsioonil ka oma nimi - tegemist on nn tensorkorrutisega või väliskorrutisega (outer product) ning R-is on selleks oma tehtemärk %o%. Nii et tuleb leida p %o% p. Selleks, et testi tulemusena korrektsed numbrid saada, tuleb esimese argumendina anda sündmuste paaride toimumise tabelist moodustatud vektor (mille saab käsuga c(maatriks)) ja tõenäosuste maatriksist samamoodi moodustatud tõenäosuste vektor kujul p=c(p_maatriks).

Märkusena tasub ehk mainida, et kui käsule chisq.test() anda ette ainult sündmusepaaride toimumise tabel maatrikskujul, siis kontrollib test samuti katsete sõltumatust, aga tabeliruutudele vastavate sündmuste tõenäosusi hinnatakse andmete põhjal (st suvandiga p= antud väärtuseid ei arvestata). Kuna vaadeldavas situatsioonis me teame, mis tõenäosusega iga sündmustepaar peaks esinema, on eelnevas toodud test tugevam.

Teema 3: Vahemike testi rakendamine sõltumatuse kontrollimiseks

Vahemike test:

Fikseerime mingi sündmuse \(A\), mille toimumise tõenäosus ühel katsel on \(p\)
Vaatleme, kui palju katseid kulub, et pärast eelmist \(A\) toimumist see toimuks jälle, moodustame jada nendest arvudest
Ooteaegade jaotuseks on \(Geom(p)\). Kirjanduses on sellest jaotusest kaks versiooni, mis erinevad vähima võimaliku väärtuse osas (kas väärtused hakkavad nullist või ühest). Meie vaatame versiooni, mis hakkab ühest, st \(P(\{X=k\})=p\cdot (1-p)^{k-1},\ k=1,2,\ldots\) ja \(P(\{X\geq \ell\})=(1-p)^{\ell-1}\)
Valime \(\ell\), mis on oluliselt väiksem kui \(n\) , võrdleme pikkustega \(1,2,...,\ell-1\) ja pikkusega vähemalt \(\ell\) vahemike arvu oodatava arvuga \(k\)-s katses (jälle Hii-ruut test)

Harjutus 4

Genereeri 10000 pseudujuhuslikku arvu lineaarse kongruentse geneaatoriga juhul \(a=5\), \(b=11\) ja \(m=2^{30}\), kasutades seemet 11. Vaatleme sündmust \(A\), mis vastab saadud katsetulemuse kuulumisele lõiku \([0,\, 0.1]\). Kontrolli vahemike testiga tulemuste vastavust katsete sõltumatuse tingimusele juhul \(\ell=11\). Mis on otsus olulisuse nivool \(\alpha=0.1\)?

Vihjed ülesande lahendamiseks:

Tekita TRUE/FALSE väärtustega vektor, mis annab iga pseudojuhusliku arvu korral teada, kas vaadeldav sündmus toimus või mitte (seda saab teha ühe käsuga, ilma tsükliteta, kasutades vektori võrdlemist arvuga). 2) Seejärel saab which käsuga leida toimumiste ajad (nt kui on vektor kas_toimus väärtustega FALSE,TRUE,FALSE,FALSE,TRUE, siis millal_toimus <- which(kas_toimus) annab tulemuseks 2,5
Kui ooteaegadele ette panna c() käsuga arv 0, siis diff() käsuga saab leida vektori ooteaegadest (eelmises näites saame which(millal_toimus) abil ooteaegadeks 2,3). Ülesande lahendamisel leitud ooteaegadest esimesed 5 peaks olema 1, 1, 1, 1, 1).
Tabeli tekitamiseks muutke seejärel kõik 11-st suuremad ajad võrseks arvuga 11,mida saab jällegi teha ühe käsuga kuju vektor[tingimus]<- uus_väärtus. NB! Selleks, et sagedustabeli tekitamisel oleks esindatud kõik vaadeldavad väärtused isegi juhul, kui vastava pikkusega ooteaegu reaalselt ei esinenud, tuleks enne tabeli tekitamist leitud ooteagade vektori faktortunnuseks käsuga kujul ajad <- factor(ajad, levels=1:ell)
Lõpuks tuleb Hii-ruut testi rakendamiseks leida ooteaegade sagedustabel ning geomeetrilise jaotuse põhjal vastavate tõenäosuste vektor. Testi rakendamisel peaks tulemuseks olema

## 
##  Chi-squared test for given probabilities
## 
## data:  tabel
## X-squared = 149.23, df = 10, p-value < 0.00000000000000022

Teema 4: Seeriate test sõltumatuse kontrollimiseks

Seeriate test

Fikseerime \(m\) vastastikku välistavat võrdvõimalikku sündmust \(A_1,\ldots, A_m\)
vaatleme, kui mitme katse järel toimub eelnevatest erinev sündmust (st vaatame, kui pikk sama sündmuse jada toimus ennem muutumist)
Võrdleme erinevate seeriate pikkuste esinemissagedusi teoreetiliste tõenäosustega, milleks on \(Geom(\frac{m-1}{m})\)

Harjutus 5

Tekita käsuga rnorm 10000 pseudojuhuslikku arvu standardsest normaaljaotusest. Vaatleme kümmet sündmust, mis vastab reaaltelje jagamisele 0.1, 0.2,…,0.9 kvantiilidega osadeks. Kontrolli vahemike testiga tulemuste vastavust sõltumatutele katsetele, vaadeldes seeriate pikkusi 1,2,3,4, ja neljast pikemaid seeriaid. Kas \(\chi^2\)-jaotuse kasutamine \(p\)-väärtuste arvutamiseks on õigustatud? Mis on tekitatud arvude korral otsus olulisuse nivool \(\alpha=0.01\)?

Vihjed lahendamiseks. Kvantiilid saab leida käsuga qnorm. Leitud kvantiile ja cut käsku kasutades saab teha kindlaks, milline sündmus igale arvule vastab. Kui muutuda cut kasu tulemus käsuga as.numeric() täisarvude vektoriks (st saame sündmuste numbrid), siis edasi saame diff käsu tulemuse ja whichkäsu abil leida kohad, kus üks sündmus läheb teiseks üle. Edasi on tegevused juba tuttavad. Kui selgub, et \(\chi^2\)-jaotuse kasutamine ei ole teststatistiku väärtuse arvutamisel õigustatud, saab R käsus chisq.test() anda ette suvandi simulate.p.value=TRUE, mille tulemusena leitakse teststatistiku \(p\)-väärtus arvutisimulatsioonide teel.

Kontrollimiseks: kui otse enne käsu rnorm kasutamist anda käsk set.seed(2023), siis esimesele kolmele x väärtusele vastavad sündmuste numbrid on 5, 2, 1 ja esimesed 16 ooteaega on 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 2, 1, 2, 1, 1