Ülesanne 2.1

Defineeri funktsioon, mis etteantud kahekomponedilisest seemnest (x₀, x₁) ja positiivsest täisarvust n lähtuvalt genereeriks n pseudojuhuslikku arvu vastavalt valemile u_i−1 = x_i / 2³⁰ , i = 2, … , n + 1,
kus arvud x_i leitakse vastavalt valemile x_i = (5 * (x_i−1)² + 171 * x_i−2 + 3) mod 2³⁰

generaator <- function(x0, x1, n) {
  u <- rep(NA, n)
  
  xim1 <- x1
  xim2 <- x0
  
  for (i in 2:(n+1)) {
    xi <- (5 * xim1^2 + 171 * xim2 + 3) %% 2^30
    u[i-1] <- xi / 2^30
    
    xim2 <- xim1
    xim1 <- xi
  }
  
  return(u)
}

Kontrolli n=10000 ja seemne (1974, 42) korral, kas sündmuse 0.2 ≤ U ≤ 0.3 toimumine on kooskõlas ühtlase jaotuse eeldusega.

n <- 10000
u <- generaator(1974, 42, n)

Fikseerime sündmuse A = {0.2 ≤ U ≤ 0.3}. Sündmuse A toimumiste arv n sõltumatu katse jooksul (N_A).

N_A = sum((u >= 0.2) * (u <= 0.3))
N_A

## [1] 1005

N_A on binoomjaotusega parameetritega n ja p, kus p on sündmuse A toimumise tõenäosus ühel katsel. Antud juhul on p lõigu [0.2; 0.3] pikkus ehk 0.1. Kvantiili definitsioonist tulenevalt peaks juhuslik suurus oma 2,5% ja 97,5% kvantiilide vahele jääma 95% tõenäosusega.

ülempiir <- qbinom(0.975, n, 0.1)
alampiir <- qbinom(0.025, n, 0.1)
c(alampiir,ülempiir)

## [1]  942 1059

(N_A >= alampiir) && (N_A <= ülempiir)

## [1] TRUE

Päris juhuslike katsete korral peaks leitud arv jääma tõenäosusega 0.05 ka piiridest välja.

Korda katset 200 korda (n = 1000, võtteks esimeseks seemneks (1974, 42) ja järgmistel kordadel viimasest generaatori väljastatud arvudele u_n ja u_n-1 vastav x_n ja x_n-1) ning leia, mitu korda jäi leitud tulemus piiridest välja.

kordi <- 200
n <- 1000

x0 <- 1974
x1 <- 42

#uued alampiir ja ülempiir (n = 1000)
ülempiir <- qbinom(0.975, n, 0.1)
alampiir <- qbinom(0.025, n, 0.1)

väljas <- 0

for(i in 1:200) {
  u <- generaator(x0, x1, n)
  N_A <- sum((u >= 0.2) * (u <= 0.3))
  
  if(N_A > ülempiir | N_A < alampiir) {
    väljas <- väljas + 1
  }
  
  x0 <- round(2^30 * u[n - 1])
  x1 <- round(2^30 * u[n])
}

väljas

## [1] 9

Kordade arv 200 katse seast, millel N_A jääb kvantiilide vahemikust välja, on ka omakorda binoomjaotusega parameetritega n = 200 (see on katsete arv) ja p = 0.05 (kvantiili definitsioonist tulenevalt peaks juhuslik suurus oma 2.5% ja 97.5% kvantiilide vahelt välja jääma 5% tõenäosusega).

Selle binoomjaotuse keskväärtus on np = 200 · 0, 05 = 10. Seega oodatavalt võiks 200-st katsest N_A jääda kvantiilide vahemikust välja 10 korda.

vahe <- abs(väljas - 10)
vahe

## [1] 1

Praegu jäi suurus N_A kvantiilide vahemikust välja 9 korral. See erineb oodatavast välja jäämiste arvust 1 võrra.

Kui suur on tõenäosus, et see kordade arv erineks oodatavast väljajäämiste arvust vähemalt niipalju, kui praguse katse tulemus näitas?

P(N_A erineb oodatavast tulemusest 10 vähemalt sama palju kui meie konkreetne tulemus 9)
= P(|N_A − 10| ⩾ vahe)
= P(10 − vahe ⩾ N_A ⩾ 10 + vahe)
= P({N_A ⩽ 10 − vahe} ∪ {N_A > 10 + vahe − 1})
= P({N_A ⩽ 10 − vahe}) + 1 − P({N_A ⩽ 9 + vahe})
= F_{N_A} (10 − vahe) + 1 − F_{N_A} (9 + vahe)
= F_{N_A} (9) + 1 − F_{N_A} (10)

Kuna N_A ∼ B(200; 0.05), siis saame F_{N_A} arvutamisel kasutada R funktsiooni pbinom (binoomjaotuse jaotusfunktsioon).

if(vahe != 0){
  pbinom(10 - vahe, kordi, 0.05) + (1 - pbinom(9 + vahe, kordi, 0.05))
} else {1}

## [1] 0.8716426

Tõenäosus, et see kordade arv erineks oodatavast väljajäämiste arvust vähemalt niipalju, kui praguse katse tulemus näitas, on ~0.87.

Ülesanne 2.2

Kirjuta Wichmann-Hilli generaator R funktsioonina, mis võtaks argumentideks arvu, mitu juhuslikku suurust on vaja teha, ja seemne (3 arvu).

WH_generaator <- function(n, xi, yi, zi) {
  tulemus <- rep(NA, n)
  for (i in 1:(n)) {
    xi <- (171 * xi) %% 30269 
    yi <- (172 * yi) %% 30307
    zi <- (170 * zi) %% 30323
    tulemus[i] <- (xi/30269 + yi/30307 + zi/30323) %% 1
  }
  return(tulemus)
}

Genereeri selle generaatoriga 1000 pseudojuhuslikku arvu, lähtudes seemnest (x0, y0, z0) = (166, 9, 32) ning kirjeldada joonisel nende arvude sattumist vahemikesse [0, 0.3), [0.3, 0.8), [0.8, 1).

Samal joonisel kujutada teoreetilised 0.05 ja 0.95 kvantiilid (mis on iga vahemiku jaoks erinevad). Kommenteeri tulemust - kas see katse tekitab mingeid kahtluseid generaatori korrektsuse osas või mitte.

n <- 1000
u <- WH_generaator(n, 166, 9, 32)

Teoreetilised kvantiilid [0, 0.3) vahemiku jaoks.

x1 <- c(0, 0.3)
alampiir1 <- qbinom(0.05, n, 0.3)
ülempiir1 <- qbinom(0.95, n, 0.3)
y1a <- c(alampiir1, alampiir1)
y1ü <- c(ülempiir1, ülempiir1)

c(alampiir1, ülempiir1)

## [1] 276 324

Teoreetilised kvantiilid [0.3, 0.8) vahemiku jaoks.

x2 <- c(0.3, 0.8)
alampiir2 <- qbinom(0.05, n, 0.8 - 0.3)
ülempiir2 <- qbinom(0.95, n, 0.8 - 0.3)
y2a <- c(alampiir2, alampiir2)
y2ü <- c(ülempiir2, ülempiir2)

c(alampiir2, ülempiir2)

## [1] 474 526

Teoreetilised kvantiilid [0.8, 1) vahemiku jaoks.

x3 <- c(0.8, 1)
alampiir3 <- qbinom(0.05, n, 1 - 0.8)
ülempiir3 <- qbinom(0.95, n, 1 - 0.8)
y3a <- c(alampiir3, alampiir3)
y3ü <- c(ülempiir3, ülempiir3)

c(alampiir3, ülempiir3)

## [1] 179 221

Histogramm ja teoreetilised 0.05 ja 0.95 kvantiilid samal joonisel.

hist(u, breaks = c(0, 0.3, 0.8, 1), ylim = c(0, 600), freq = TRUE)

lines(x1, y1a, col = "red", lty = 3)
lines(x1, y1ü, col = "red", lty = 3)

lines(x2, y2a, col = "red", lty = 3)
lines(x2, y2ü, col = "red", lty = 3)

lines(x3, y3a, col = "red", lty = 3)
lines(x3, y3ü, col = "red", lty = 3)

See katse tekitab kahtluseid generaatori korrektsuse osas, kuna graafikul näeme et esimeses vahemikus on vähem arve, kui näitab 0.05 kvantiil.

Aga me ei saa öelda, kas see geraator on hea või mitte, sest juhuslik suurus peab oma 0.05 ja 0.95 kvantiilide vahelt välja jääma 10% tõenäosusega. Järelikult on vaja rohkem katseid.