Homework 4

Aleksei Parm

2022-10-03

Etteantud jaotusega pseudojuhuslike arvude genereerimine jaotusfunktsiooni pööramise meetodil

Tekita vaikeparameetritega runif generaatoriga seemnest 199 lähtuvalt 100000 pseudojuhuslikku arvu ning kasuta neid selleks, et saada jaotusfunktsioonile

vastavad pseudojuhuslikud arvud.

NB! Kindlasti veendu, et leitud pöördfunktsiooni korral kehtib võrdus F(F−1(u)) = u ∀u ∈ (0, 1) (seda võib teha graafiku joonistamise teel).


Tekitame pseudojuhuslikud arvud.

n <- 100000
set.seed(199)
u <- runif(n)


Defineerime jaotusfunktsiooni ja selle pöördfunktsiooni.

Pöördfunktsiooni leidmiseks tuleb lahendada võrrand F(x) = u muutuja x suhtes.

Fjaotus <- function(x){
  return(ifelse(x < 0, 0, ifelse(x < 3, x/6, (x-2)^5 / (1+(x-2)^5))))
}

Finv <- function(u){
  return(ifelse(u < 0, 0, ifelse(u < 1/2, 6*u, (u/(1-u))^(1/5) + 2)))
}

curve(Fjaotus(Finv(x)), 0.01, 0.99)

Graafikul näeme, et leitud pöördfunktsiooni korral kehtib võrdus F(F−1(u)) = u ∀u ∈ (0, 1)

x = Finv(u)

Testid

Kontrolli saadud arvude vastavust toodud jaotusfunktsioonile Cramer-von Mises testiga ning sõltumatust seeriate testiga (vaadeldes seeriaid pikkusega 1,2,3,4 ja neljast suuremaid), defineerides kvantiilide abil 9 võrdse tõenäosusega sündmuse täissüsteemi.

Cramer-von Mises

#install.packages("goftest")
goftest::cvm.test(x, Fjaotus)
## 
##  Cramer-von Mises test of goodness-of-fit
##  Null hypothesis: distribution 'Fjaotus'
##  Parameters assumed to be fixed
## 
## data:  x
## omega2 = 0.04581, p-value = 0.9012

p-väärtus on väga suur (p-value = 0.9012 >> 0.1).

Jääme nullhüpoteesi juurde (arvud vastavad jaotusfunktsiooniga F juhuslikele suurustele)


Seeriate test (sõltumatuse kontrollimiseks)

Leiame jaotuspunktideks kasutatavad reaalarvud jaotusfunktsiooni pöördfunktsiooni kasutades:

cp <- c(-Inf, Finv((1:9)/10),Inf)
cp
##  [1]     -Inf 0.600000 1.200000 1.800000 2.400000 3.000000 3.084472 3.184664
##  [9] 3.319508 3.551846      Inf


Kasutame cut käsku, et teha kindlaks, mitmendasse kvantiilide poolt tekitatud osalõiku iga arv kuulub:

sündmused <- as.numeric(cut(x, breaks=cp))
sündmused[1:10]
##  [1] 1 7 8 3 1 3 7 7 3 4


Uus seeria algab, kui vaadeldaval kohal olev sündmus erineb eelmisest ehk toimunud sündmuste numbrite vahe erineb nullist. Võime mõelda nii, et uus seeria algab ka esimesest vaatlusest, siis enne vahede leidmist kirjutame ette arvu, mis jadas ei saa esineda

algused = which(diff(c(-1, sündmused)) != 0)
algused[1:10]
##  [1]  1  2  3  4  5  6  7  9 10 11


Seeriate pikkused on algusaegade vahed.

seeriad = diff(c(algused))
seeriad[1:10]
##  [1] 1 1 1 1 1 1 2 1 1 1


Edasi käib kõik samamoodi nagu varem, ainult et sündmuse(uue seeria algamise) tõenäosus on 0.9 ja kokku paneme kõik neljast pikemad seeriad.

Kuna neljast pikemate seeriate saamise tõenäosus on ainult 0.0001, siis peab Hii-ruut testi kasutamisel olema ettevaatlik. Statistiku väärtusele vastava p-väärtuse arvutamisel ei ole χ2 jaotuse kasutamine õigustatud, seega tuleks p väärtus leida simulatsioonide teel.

ell = 5
seeriad[seeriad > ell] = ell
tabel = table(factor(seeriad, levels = 1:ell))
pA = 0.9
p <- c(pA * (1-pA)^(0:(ell-2)), (1-pA)^(ell-1))
chisq.test(tabel, p=p, simulate.p.value = TRUE)
## 
##  Chi-squared test for given probabilities with simulated p-value (based
##  on 2000 replicates)
## 
## data:  tabel
## X-squared = 1.3204, df = NA, p-value = 0.8576

p-väärtus >> 0.1, jääme nullhüpoteesi juurde (genereeritud pseudojuhuslikud arvud käituvad nagu sõltumatute katsete tulemused)