[vastaus aiempaan viestiin]
| Kirjoittaja: | Seppo Mustonen |
|---|---|
| Sähköposti: | - |
| Päiväys: | 8.11.2004 14:20 |
Katsellessani verkosta, mitä viime aikoina on tapahtunut satunnaisluku-
tutkimuksen alueella, huomasin erään uuden generaattorin saaneen
melkoista suosiota ja päätin liittää sen nykyiseen Survoon.
Kyseessä on kahden japanilaisen (Takuji Nishimura ja Makoto Matsumoto)
vuonna 1998 esittämä (ja vuonna 2002 parantama) "Mersenne Twister",
jolla on erittäin suuri jakso M=2^19937-1 eli M on lähes 10^6002.
Generaattorin perustana oleva teoria on kuvattu tekijöiden artikkelissa
"Mersenne Twister: A 623-dimensionally equidistributed uniform
pseudorandom number generator" ja luettavissa osoitteesta
http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/ARTICLES/mt.pdf
(aika rankkaa teoriaa).
Keskeinen parametri on alkuluku M (vuonna 1971 löydetty Mersennen
luku). Mersennen luvuista tarkemmin:
http://www.utm.edu/research/primes/mersenne/
Uusi generaattori on saatavilla funktiona mrand() SURVO MM:n
operaatioissa VAR, TRANSFORM BY #UNIFORM, MAT #TRANSFORM ja
MNSIMUL versiosta 2.18 lähtien.
* * *
Hyvien pseudosatunnaislukujen tekeminen ei ole helppoa.
Osoitan pienellä esimerkillä, miten omituisia tilanteita voi syntyä,
vaikka satunnaislukuja tuotetaan periaatteessa kelvollisella
generaattorilla.
Vanhemmassa kirjallisuudessa yleisimmin on esitelty ns. sekakongruenssi-
menetelmää, jossa peräkkäiset satunnaisluvut U(0),U(1),U(2),...,U(n),...
muodostetaan kaavalla
U(n)=mod(a*U(n-1)+c,m) , U(0)=siemenluku.
Tässä mod(N,M) tarkoittaa jakojäännöstä, kun kokonaisluku N jaetaan
konaisluvulla M, esim. mod(26,2)=0 ja mod(26,4)=2.
Alunperin kaikki tapahtuu kokonaisluvuilla (jolloin saadaan samoja
tuloksia riippumatta koneen aritmetiikkaominaisuuksista) ja "tasaista
välin [0,1) jakaumaa" noudattaviksi tarjotaan silloin lukuja U(n)/m.
..........
Kun a=501, c=9^5, m=10^5 ja U(0)=0, niin Survon laskentakaavio
U(N)|=if(N=0)then(0)else(mod(a*U(N-1)+c,m))
tuottaa seuraavia lukuja
U(0)=0
U(1)=59049
U(2)=42598
U(3)=647
U(4)=83196
U(5)=40245
U(6)=21794
U(7)=77843
U(8)=58392
U(9)=13441
Tiedetään, että tällainen generaattori antaa kaikki kokonaisluvut
0,1,2,...,m-1 jossain järjestyksessä (eli jakson pituus on m), mikäli
1) luvuilla c ja m ei ole yhteisiä tekijöitä,
2) mod(a,p)=1 kaikilla luvun m alkulukutekijöillä p,
3) mod(a,4)=1.
Nämä ehdot eivät kuitenkaan takaa välttämättä mitään satunnaisuuden
tapaista. Esim. valinta a=1, c=1 täyttää ehdot 1-3, mutta synnyttää
nollasta lähtien peräkkäin luvut 0,1,2,...,m-1,0,1,2,...,m-1,0,1,2,...
oli m mikä tahansa.
Edellä kuvattu valinta a=501, c=9^5, m=10^5 ja U(0)=0 saattaa vaikuttaa
paremmalta, mutta kohta nähdään, etteivät luvut tule silloinkaan olemaan
kelvollisia.
Havainnollisuuden vuoksi osoitan ensin Survon avulla,
että po. tapauksessa syntyvät kaikki luvut 0,1,2,...,m-1
jossain järjestyksessä, ennenkuin samat luvut alkavat toistua.
Tämähän on etukäteenkin selvää, sillä "jaottomuusehdot" 1-3
ovat voimassa.
..........
Luodaan kahden muuttujan U,X datatiedosto JONO:
FILE CREATE JONO,16,2
FIELDS:
1 N 8 U (#####)
2 N 8 X
END
Varataan tilaa m+1 havainnolle:
FILE INIT JONO,100001
Tehdään luvut U(0),U(1),...,U(m-1),U(m):
VAR U TO JONO
U=if(ORDER=1)then(0)else(U2)
U2=mod(a*U[-1]+c,m) a=501 c=9^5 m=10^5
..........
Otetaan esiin 10 ensimmäistä U-lukua:
FILE LOAD -JONO,CUR+1 / IND=ORDER,1,10 VARS=U
0
59049
42598
647
83196
40245
21794
77843
58392
13441
Ne ovat samoja kuin aikaisemmat laskentakaaviolla tehdyt.
..........
Viimeisen 100001. luvun tulisi olla sama kuin ensimmäinen U(0)=0,
koska jakso alkaa uudelleen. Näin on todella asian laita, sillä
FILE LOAD -JONO,CUR+1 / IND=ORDER,100001,100001 VARS=U
0
..........
Tämä ei vielä riitä osoittamaan, että ensimmäiset 100000 lukua
ovat luvut 0,1,2,...,99999 jossain järjestyksessä.
Siksi asetetaan luvut U suuruusjärjestyksessä tiedostoon JONOS:
FILE SORT JONO BY U TO JONOS / IND=ORDER,1,100000
..........
Alku näyttää hyvältä:
FILE LOAD -JONOS,CUR+1 / IND=ORDER,1,10 VARS=U
0
1
2
3
4
5
6
7
8
9
..........
100000 luvun tarkistaminen silmämääräisesti on turhauttavaa. Siksi
kannattaa käyttää laskennallisia keinoja.
Jos kaikki 100000 lukua ovat todella mukana ja suuruusjärjestyksessä,
peräkkäisten lukujen erotus on aina 1.
Lasketaan nuo erotukset:
VAR X=if(ORDER=1)then(1)else(U-U[-1]) TO JONOS
..........
Kaikkien X-arvojen tulisi olla ykkösiä. Tarkistetaan:
STAT JONOS,CUR+1 / VARS=X
Basic statistics: JONOS N=100000
Variable: X ~if(ORDER=1)then(1)else(U-U[-1])
Constant= 1
X on siis vakio 1. Saatiin odotettu tulos eli generaattorin jakso
on todella m=10^5.
..........
Tutkitaan nyt hieman lisää generaattorin ominaisuuksia.
Munnetaan luvut U välin [0,1) luvuiksi X:
VAR X=U/100000 TO JONO / IND=ORDER,1,100000
..........
On useita tapoja osoittaa, että X-luvut eivät ole riittävän riippumatto-
mia toisistaan. Tarkastellaan peräkkäisten lukujen muodostamia pareja
X(n),X(n-1), n=2,...,m.
Määritellään U uudelleen viivästettynä X(n)-arvona X(n-1):
VAR U=X[-1] TO JONO / IND=ORDER,2,100000
..........
Piirretään hajontakuva 4000 ensimmäisestä pisteparista X(n),X(n-1):
GPLOT JONO,X,U / IND=ORDER,2,4001
Pisteiden tulisi sijoittua suhteellisen tasaisena mattona yksikkö-
neliöön. Näin ei tässä tapauksessa käy, vaan matossa esiintyy hyvin
systemaattista kuviointia ja siihen jää paljon valkoisia läiskiä.
Tämä on tyypillinen sekakongruenssimenetelmän heikkous.
Kun jakso m on riittävän suuri ja vakiot valitaan taitavasti, nämä
puutteet eivät tule yhtä selvästi esille.
..........
Juuri nyt tarkastellulla generaattorilla on myös seuraava todella
kummallinen ominaisuus.
Talletetaan 400 ensimmäistä X-arvoa vektoriksi, joka muunnetaan
20*20-matriisiksi A:
MAT SAVE DATA JONO TO X / VARS=X IND=ORDER,1,400
MAT A=VEC(X,20) / *A~VEC(X) 20*20
Poimitaan näytteeksi "satunnaismatriisin" A 8 ensimmäistä saraketta:
MAT LOAD A(1:20,1:8),1.12345,CUR+1
MATRIX A
VEC(X)
/// 1 2 3 4 5 6 7 8
1 0.00000 0.35980 0.71960 0.07940 0.43920 0.79900 0.15880 0.51860
2 0.59049 0.85029 0.11009 0.36989 0.62969 0.88949 0.14929 0.40909
3 0.42598 0.58578 0.74558 0.90538 0.06518 0.22498 0.38478 0.54458
4 0.00647 0.06627 0.12607 0.18587 0.24567 0.30547 0.36527 0.42507
5 0.83196 0.79176 0.75156 0.71136 0.67116 0.63096 0.59076 0.55056
6 0.40245 0.26225 0.12205 0.98185 0.84165 0.70145 0.56125 0.42105
7 0.21794 0.97774 0.73754 0.49734 0.25714 0.01694 0.77674 0.53654
8 0.77843 0.43823 0.09803 0.75783 0.41763 0.07743 0.73723 0.39703
9 0.58392 0.14372 0.70352 0.26332 0.82312 0.38292 0.94272 0.50252
10 0.13441 0.59421 0.05401 0.51381 0.97361 0.43341 0.89321 0.35301
11 0.92990 0.28970 0.64950 0.00930 0.36910 0.72890 0.08870 0.44850
12 0.47039 0.73019 0.98999 0.24979 0.50959 0.76939 0.02919 0.28899
13 0.25588 0.41568 0.57548 0.73528 0.89508 0.05488 0.21468 0.37448
14 0.78637 0.84617 0.90597 0.96577 0.02557 0.08537 0.14517 0.20497
15 0.56186 0.52166 0.48146 0.44126 0.40106 0.36086 0.32066 0.28046
16 0.08235 0.94215 0.80195 0.66175 0.52155 0.38135 0.24115 0.10095
17 0.84784 0.60764 0.36744 0.12724 0.88704 0.64684 0.40664 0.16644
18 0.35833 0.01813 0.67793 0.33773 0.99753 0.65733 0.31713 0.97693
19 0.11382 0.67362 0.23342 0.79322 0.35302 0.91282 0.47262 0.03242
20 0.61431 0.07411 0.53391 0.99371 0.45351 0.91331 0.37311 0.83291
Matriisi A vaikuttaa hyvin umpimähkäiseltä (myös tässä näkymättömien
alkioiden osalta).
Lasketaan matriisin A käänteismatriisi B:
MAT B=INV(A,det) / *B~INV(VEC(X)) det=-0.00742896 20*20
Matriisin A determinantti on aika pieni, mutta ei silti aiheuta
numeerisia ongelmia matriisin kääntämisessä.
Tarkistetaan käänteismatriisin kelvollisuus laskemalla C=A*B, jonka
pitäisi olla riittävän tarkkaan yksikkömatriisi:
MAT C=A*B / *C~VEC(X)*INV(VEC(X)) 20*20
MAT E=IDN(20,20)-C / Yksikkömatriisin ja C:n erotusmatriisi
MAT S=SUM(SUM(E,2)') / Erotusmatriisin alkioiden neliösumma
MAT_S(1)=1.2956007578092e-026
Erotusmatriisin alkioiden neliösumma osoittaa, että matriisin A
käänteismatriisi on B on suhteellisen tarkka.
Käänteismatriisi B näyttää kuitenkin tällaiselta:
LOADM B,(C5),CUR+1
INV(VEC(X))
1 2 3 4 5 6 7 8 9 10
1 6.0 0.00 0.00 -0.00 0.49 0.00 0.00 1.00 -0.00 -0.00
2 -25.5 1.00 -1.00 0.00 1.05 -1.00 -2.00 -3.00 0.00 0.00
3 31.4 -1.00 1.00 -0.00 -0.56 1.00 3.00 3.00 -0.00 -0.00
4 -6.4 -0.00 -0.00 0.00 0.51 -0.00 -1.00 -1.00 0.00 0.00
5 12.2 0.00 0.00 -0.00 -0.02 1.00 1.00 1.00 -0.00 -0.00
6 12.5 0.00 0.00 -0.00 -0.02 0.00 1.00 1.00 -0.00 -0.00
7 -6.1 -0.00 -0.00 0.00 0.51 -1.00 -1.00 -1.00 1.00 1.00
8 -19.8 -0.00 -0.00 0.00 0.54 -0.00 -1.00 -1.00 -1.00 -1.00
9 -24.5 -0.00 -0.00 0.00 1.05 -1.00 -3.00 -3.00 0.00 0.00
10 17.3 0.00 1.00 -1.00 -2.54 1.00 3.00 3.00 -0.00 -0.00
11 19.0 0.00 0.00 -0.00 -0.54 0.00 2.00 2.00 -0.00 -0.00
12 -26.7 -0.00 -1.00 0.00 0.05 -0.00 -2.00 -2.00 0.00 0.00
13 43.3 -1.00 1.00 1.00 -0.59 2.00 4.00 5.00 -1.00 -0.00
14 -36.7 1.00 -1.00 0.00 1.07 -2.00 -4.00 -5.00 1.00 0.00
15 -7.4 -0.00 -0.00 0.00 -0.49 -0.00 -0.00 -1.00 0.00 0.00
16 -6.5 -0.00 -0.00 0.00 0.51 -0.00 -1.00 -0.00 0.00 -1.00
17 19.1 0.00 0.00 1.00 0.46 0.00 1.00 1.00 -0.00 1.00
18 -32.1 -0.00 -1.00 0.00 0.56 -1.00 -3.00 -3.00 0.00 0.00
19 6.8 0.00 0.00 -0.00 -0.51 0.00 1.00 0.00 -0.00 -0.00
20 18.6 0.00 1.00 -1.00 -1.54 1.00 2.00 3.00 -0.00 -0.00
11 12 13 14 15 16 17 18 19 20
1 1.00 0.00 6.0 -0.00 0.00 -6.0 -0.00 -6.0 -0.00 -0.00
2 -3.00 -2.00 -26.5 2.00 1.00 25.5 0.00 25.5 -1.00 -1.00
3 3.00 3.00 32.4 -2.00 -1.00 -31.4 -0.00 -31.4 1.00 1.00
4 -1.00 -1.00 -6.4 1.00 -0.00 6.4 0.00 6.4 0.00 0.00
5 1.00 1.00 13.2 -1.00 0.00 -12.2 -0.00 -12.2 -0.00 -0.00
6 1.00 1.00 12.5 -1.00 0.00 -12.5 -0.00 -12.5 1.00 1.00
7 -1.00 -1.00 -7.1 1.00 -0.00 6.1 0.00 6.1 0.00 0.00
8 -2.00 -1.00 -20.8 0.00 1.00 19.8 0.00 20.8 -1.00 -1.00
9 -3.00 -3.00 -26.5 2.00 1.00 25.5 1.00 25.5 -1.00 -1.00
10 4.00 3.00 20.3 -3.00 -2.00 -18.3 -1.00 -19.3 1.00 1.00
11 2.00 1.00 20.0 -1.00 -1.00 -19.0 -0.00 -19.0 1.00 1.00
12 -1.00 -1.00 -26.7 1.00 -0.00 26.7 -1.00 25.7 -1.00 -1.00
13 4.00 4.00 44.3 -3.00 -1.00 -43.3 -0.00 -43.3 1.00 1.00
14 -5.00 -4.00 -38.7 4.00 1.00 36.7 1.00 37.7 -1.00 -1.00
15 -0.00 -0.00 -6.4 0.00 -0.00 6.4 0.00 6.4 0.00 0.00
16 -1.00 -1.00 -7.5 0.00 1.00 7.5 0.00 7.5 0.00 -1.00
17 1.00 1.00 19.1 -0.00 0.00 -19.1 -0.00 -19.1 -0.00 1.00
18 -3.00 -2.00 -33.1 2.00 1.00 32.1 0.00 32.1 -1.00 -1.00
19 0.00 0.00 6.8 -0.00 0.00 -6.8 1.00 -6.8 -0.00 1.00
20 3.00 2.00 19.6 -2.00 -1.00 -18.6 -1.00 -18.6 1.00 -0.00
Matriisin B sarakkeista vain 5 (20:stä) on umpimähkäisen oloisia.
Loput 15 saraketta ovat muodostuneet pienistä kokonaisluvuista!
Koska epäilin, että joillakin pyöristyksillä olisi osuutta tuloksiin,
käänsin kokonaislukumatriisin 100000*A tarkasti Maple-ohjelmalla ja
sain yhtäpitävät tulokset. Itse asiassa likimääräiset kokonaisluvut
ovat tarkkoja. Esim. käänteismatriisin inv(A) ensimmäinen vaakarivi
on tarkasti laskettuna
[355407 98
[------ , 0 , 0 , 0 , --- , 0 , 0 , 1 , 0 , 0 , 1 , 0 ,
[58960 201
355407 -355407 -355407 ]
------ , 0 , 0 , ------- , 0 , ------- , 0 , 0]
58960 58960 58960 ]
ja mm. vastaavuus numeerisesti laskettuun B-matriisiin on hyvä:
MAT_B(1,1)=6.027934192673
355407/58960=6.027934192673
MAT_B(1,5)=0.48756218905473
98/201=0.48756218905473
Tämä on aika omituinen juttu. On kyllä tyypillistä, että
käännettäessä matriisi, jonka alkiot ovat yksinkertaisia kokonaislukuja
saadaan tulos, jossa alkiot ovat kaikkea muuta kuin kokonaislukuja.
Enpä olisi odottanut kohtaavani tilannetta, jossa tämä ominaisuus olisi
kääntynyt niin, että "satunnaismatriisin" käänteismatriisi on
rakenteeltaan paljon yksinkertaisempi kuin alkuperäinen!
Löysin tämän omituisuuden jo vuonna 1976 ja se toimii edelleen
varoittavana esimerkkinä. Kannattaa olla todella tarkkana simulointi-
kokeita tehdessään. Millä tahansa pseudosatunnailukugeneraattorilla
voi jossain tilanteessa olla heikot hetkensä.
On paikallaan lopuksi huomauttaa siitä, että hyvälläkin generaattorilla
luodussa satunnaismatriisissa on aina jonkinlaista rakenteellista
systematiikkaa, sillä käänteismatriisin alkiot eivät ole toisistaan
riippumattomia.
Tätä piirrettä havainnollistan kääntämällä 200*200-satunnaismatriisin
ja piirtämällä sen käänteismatriisin matriisidiagrammana:
..........
p=200
MAT A=ZER(p,p)
MAT #TRANSFORM A BY mrand(8112004)
MAT B=INV(A,det)
GPLOT B.MAT / TYPE=MATRIX NORM=T MODE=1000,1000
ROWLABELS= COLUMNLABELS=
Käänteismatriisin kuva on hieman räsymaton kaltainen pysty- ja vaaka-
suorine raitoineen. Tätä kannattaa verrata kuvaan alkuperäisestä
satunnaismatriisista
GPLOT A.MAT
joka on tasaista mössöä, kuten pitääkin.
-Seppo
| Vastaukset: |
|---|
Survo-keskustelupalstan (2001-2013) viestit arkistoitiin aika ajoin sukrolla, joka automaattisesti rakensi viesteistä (yli 1600 kpl) HTML-muotoisen sivukokonaisuuden. Vuoden 2013 alusta Survo-keskustelua on jatkettu entistäkin aktiivisemmin osoitteessa forum.survo.fi. Tervetuloa mukaan!