Statistikk og Simulering

Prosjekt 1

Veke 3. Introduksjon til slumptalsgeneratorar (labøving)

5.3. Veke 3. Introduksjon til slumptalsgeneratorar (labøving)

  5.3.1 Fyrste test i Matlab
  5.3.2 Perioden
  5.3.3 Uniform fordeling
  5.3.4 Frøet i Matlabs slumptalsgenerator
  5.3.5 Statistiske testar (Ekstra)
  5.3.6 Avslutning

Den mest populære metoden for å generera slumptal på ein datamaskin er lineære kongruensgeneratorar. Dei har mange lytar, men dei har vore velkjende over lang tid, dei er enkle å implementera, og dei er effektive. Der er ein del alternativ, og mange er betre sjølv om dei heller ikkje er lytefrie. Me tek for oss den lineære kongruensgeneratoren for å illustrera nokre generelle fenomen og utfordringar.

Definisjon 14 Ein lineær rekurrens er gjeve ved formelen

s_{i} = a \cdot s_{i - 1} + c mod m,

for konstante heiltal $a$ , $c$ og $m$ . Dersom me startar med eit frø (seed) $s_{0}$ , kan me bruka rekkurensen til å generera ei fylgje av slumptal $x_{1}, x_{2}, x_{3}, \dots$ . Denne slumptalsgeneratoren er kjend som Lehmers algoritme og som en lineær kongruensgenerator.

Oppgåve 5.15 (Refleksjon) Sjå på definisjon over. Hugsar du kva $mod$ -operatoren tyder? Sjå tilbake på pensum frå Diskret Matematikk om det trengst.

5.3.1. Fyrste test i Matlab

Me kan implementera kongruensgeneratoren i Matlab som fylgjer:

1function [x] = rng25(n,x0) 
2 
3a = int32(6) ;     % Multiplicative constant 
4c = int32(7) ;     % Additive constant 
5m = int32(25) ;    % Modulus 
6persistent s       % State 
7 
8if nargin > 1, 
9   % If x0 is given, we use it as the seed 
10   s = int32(x0) ; 
11elseif isempty(s), 
12   % If no state is set, we use a hardcoded one. 
13   s = int32(11) 
14end 
15 
16% Start with an empty array 
17x = [] ; 
18 
19% Iterate n times to generate an array of random numbers 
20for i=1:n, 
21   s = mod(a*s + c, m ) ; 
22   x = [ x s ] ; 
23end

Legg merke til eit par finessar i Matlab.

1.: Kodeordet persistent let oss definera lokale variablar som funksjonen hugsar frå kall til kall. Dette svarer til static i C.
2.: Me reknar med 32-bits heiltal. Sidan Matlab elles bruker flyttal, må alt konverterast med int32()-funksjonen.
3.: Me kan sjekka om me har fått alle argumenta, ved å sjekka nargin som gjev talet på argument.

Oppgåve 5.16 Last ned fila rng25.m. Dette er den same funksjonen som er gjengjeve over. Test funksjonen på kommandolina og genererer 25 tilfeldige tal. Vel din eigen verdi for frøet $s$ .

1rng25(25,s)

Drøft med sidemannen: ser tala tilfeldige ut?

Oppgåve 5.17 Prøv funksjonen over igjen, og generer 100 tilfeldige tal. Drøft med sidemannen: ser tala tilfeldige ut?

Oppgåve 5.18 Kopier fila rng25.m og lag ein ny funksjon der du endrar den multiplikative parameteren til $a = 9$ . Kva skjer? Kva fortel dette oss om desse slumptalsgeneratorane?

Oppgåve 5.19 Lag ein tredje variant, med $a = 10$ . Kva skjer denne gongen?

5.3.2. Perioden

Definisjon 15 Perioden i ei fylgje $[x_{1}, x_{2}, \dots]$ er det minste talet $m > 0$ slik at $x_{i} = x_{i + m}$ for alle $i > k$ for ein eller annan $k$ .

Ein fylgje som har periode $m$ repeterer seg sjølv for kvart $m$ te element. Merk at der kan vera ein serie med meir enn $m$ distinkte tal i starten av fylgja, men frå eit eller anna punkt $k$ kan me sjå periodiske syklar med fast lengd.

Oppgåve 5.20 Sjå på fylgjene som du genererte med dei tre kongruensgeneratorane over (for $a = 6, 9, 10$ ). Kva periode har fylgjene? Er periode uavhengig av frøet $s$ ?

Oppgåve 5.21 Kva er den største perioden du kan tenkja deg for ein kongruensgenerator med parameter $a$ , $c$ og $m$ .

Sats 3 Den lineære kongrunsgeneratoren $s_{i} = a s_{i - 1} + c mod m$ har periode $m$ dersom, og berre dersom,

1.: $c$ og $m$ er relativt prim (mao. $hcf (c, m) = 1$ )
2.: $p$ deler $b = a - 1$ for alle primtal $p$ som deler $m$
3.: $4$ deler $b = a - 1$ dersom $4$ som deler $m$

5.3.3. Uniform fordeling

Ein god slumptalsgenerator skal gje slumptal som er uavhengige og uniformt fordelte. Dersom me ser på ein heil periode frå ein lineære kongruensgenerator, er dette aldri tilfredsstilt. Generatoren kan ikkje gjenta tal innanfor perioden, so når periode går mot slutten veit me at neste tal må vera eit av dei som me ikkje har sett før. Dersom me ser på ein liten del av perioden (og parametra $a$ , $c$ og $m$ er velvalde), kan me derimot få ein fordeling som er tilfredsstillande nær uniform og uavhengig.

Oppgåve 5.22 Lag to nye matlabfunksjonar rng1.m med $a = 219$ , $c = 7$ og $m = 32749$ og rng2.m med same $a$ og $c$ men $m = 2^{15} = 32768$ . Merk at den fyrste $m$ -verdien er eit primtal.

Me skal testa dei to slumptalsgeneratorane, og vurdera om fordelinga er uniform. Dette kan me ikkje gjera når utfallsrommet er større enn utvalet. Difor kan me ikkje sjå eit slumptal $x$ direkte for generatoren, men dersom $x$ er uniformt fordelt, so vil ogso $x mod q$ vera uniformt fordelt.

Oppgåve 5.23 Bruk rng1.m og generer ein vektor $x$ med 1000 tilfeldige tal. Lag so ein ny vektor $y$ som er $x$ redusert elementvis modulo 16. Lag eit histogram av $y$ . Dvs.

1x = rng1(25,s) 
2y = mod(x,16) 
3histogram(y, ’BinMethod’,’integers’)

Ser det ut som om tala er uniformt fordelte?

Oppgåve 5.24 Gjenta øvinga over for rng2.m. Er dette uniformt fordelt?

Testane over ser på dei fire minst signifikante bitsa i slumptala. Me kan òg sjå på dei fire mest signifikante sifra.

Oppgåve 5.25 Bruk rng1.m og generer ein vektor $x$ med 1000 tilfeldige tal. Lag so ein ny vektor $y$ som er $x$ med dei fire mest signifikante sifra, ved å dela på $2^{11}$ . Lag eit histogram av $y$ . Dvs.

1x = rng1(25,s) 
2y = x/2^11 
3histogram(y, ’BinMethod’,’integers’)

Ser det ut som om tala er uniformt fordelte?

Gjenta øvinga med rng2.m. Kva ser du?

5.3.4. Frøet i Matlabs slumptalsgenerator

Oppgåve 5.26 Finn fram cointrial-funksjonen din frå forrige veke. Simuler åtte kast med tre myntar per kast. Gjenta denne simuleringa tre gongar. Får du same resultat kvar gong?

Oppgåve 5.27 Gjennta forrige oppgåve, men køyr kommandoen rng(243) før kvart kall til cointrial. Kva skjer no? Kvifor?

I stokastisk simulering er det mogleg å setja frøet manuelt før simuleringa. Ved å bruka same frø, kan ein gjenta eksakt same simulering determistisk. Somme tider er det nyttig, men det legg òg eit stort ansvar på brukaren for å sjå til at frøet er fornuftig (tilfeldig) valgt.

5.3.5. Statistiske testar (Ekstra)

(Dette avsnittet er ein forsmak på noko som me kjem tilbake til seinare i semesteret. Dersom du ikkje rekk over det, so er det ikkje noko problem.)

Øving 5.22 tok utgangspunkt i ein hypotese, at slumptala frå rng1.m er uniformt fordelte. Dersom denne hypotesen er sann, so er det òg sant at slumptala modulo 16 er sann. Når me ser på datamaterialet, i form av eit histogram, kan me vurdera om denne hypotesen er rimeleg eller ikkje.

Dersom histogrammet ikkje ser ut som ei uniform fordeling, er det rimeleg å tru at det ikkje er generert av ei uniform fordeling, og hypotesen er usann. Me går då ut frå at slumptalsgeneratoren er dårleg.

Dersom histogrammet ser ut som ei uniform fordeling, er det rimeleg å gå ut frå at hypotesen er sann, og me generatoren har i alle fall ikkje denne dårlege eigenskapen.

Slik visuell vurderinga av datamaterialet vert mykje synsing. For å gjera ein objektiv vurdering ynskjer me enkle kvantitative svar.

Lat oss ta utgangspunkt i histogrammet igjen. Lat $[y_{1}, y_{2}, \dots, y_{n}]$ vera utvalet vårt, dvs. ei fylgje av slumptal modulo 16.

1.: Lat $F_{y}$ vera frekvensen av verdien $y$ i utvalet.
Dvs. $F_{y}$ , for $0 \leq y \leq 15$ er talet på gongar $y$ førekjem i utvalet. Histogrammet plottar $F_{y}$ for kvar $y$ .
2.: Lat $E_{y}$ vera forventingsverdien til $F_{y}$ , dersom hypotesen vår er sann.
Hypotesa seier uniform fordelinga, og då er $E_{y} = n ∕ 16$ der $n$ er storleiken på utvalet.
3.: Me reknar ut den stokastiske variabelen $G = \sum_{y = 0}^{15} \frac{{(F_{y} - E_{y})}^{2}}{E_{y}} .$

Variabelen $G$ er eit standardverkty for å samanlikna ein empirisk fordeing ( $F_{u}$ ) med ein hypotetisk fordeling (uniform i dette tilfellet). Det er lettare å sjå avvik i ein enkelt skalarvariabel $G$ , enn å sjå på heile histogrammet med seksten forskjellige frekvensar.

Oppgåve 5.28 Sjå på uttrykket for $G$ . Kva verdiar kan $G$ ta? Korleis ser histogrammet ut når $G$ tek minste mogleg verdi?

Oppgåve 5.29 Kva verdiar ventar du at $G$ har når hypotesen om uniform fordeling held? Kva når ho ikkje held?

Oppgåve 5.30 Bruk rng1.m som i førre avsnitt og lag eit utval på $n = 1000$ tilfeldige tal modulo 16. Finn frekvensane $F_{y}$ vha. funksjonen

1f = histcounts(y,’BinMethod’,’integers’)

Rekna ut

G

som forklart over. Kva verdi får du?

Gjer det same for rng2.m.

Oppgåve 5.31 Variabel $G$ er stokastisk med $χ^{2}$ -fordeling med 15 fridomsgradar. Plott sannsynsfordelinga

1fplot( @(x)chi2pdf(x,15), [0 40] )

Det merkelege uttrykket @(x)chi2pdf(x,15) er eit lambdauttrykk og lagar ein ny funksjon med ein parameter

x

vha. den eksisterande funksjonen som har 2.

Samanlikna dine observasjonar av $G$ frå forrige oppgåve med sannsynsfordelinga. Synest du observasjonane dine ser sannsynlege ut dersom slumptala er uniformt fordelte?

5.3.6. Avslutning

I denne økta har me berre skrapt i overflata på eit stort og viktig tema. Dei generatorane som me har sett på er ganske primitive, og det er tydeleg at dei ikkje er perfekte. Det er viktig å hugsa på at dei same problema går igjen i meir sofistikerte slumptalsgeneratorar. Dei er aldri perfekte, og det vesentlege spørsmålet vil alltid vera om dei er gode nok for eit bestemt bruksområde.

Slumptalsgeneratorar er eit kontroversielt tema. Der er delte meiningar i literaturen om kva som er viktig og om kva som er godt nok.

Det er viktig ikkje å vera naiv i valet av slumptalsgeneratorar. Mange standardbibliotek har implementert generatorar med dårlege statistiske eigenskapar, sjølv om det kanskje ikkje er like ille som det var for ein generasjon sidan. På ein del aktuelle område er dei statistiske eigenskapane kritiske. Tenk deg nettpoker der korta ikkje er uniformt fordelte? Eller eit nettbanksystem der bandittar veit at nokre nyklar er meir sannsynlege enn andre? Eller ein vitskapleg simulering der sannsynsfordelinga er ei heilt anna enn det som forskaren føreset?

Me avslutta økta med eit døme på statistisk hypotesetest. Det er eit sentralt tema som me går vidare med neste økt. Hypotesetestar vert brukt på mange andre problem enn kvaliteten til slumptalsgeneratorar.