Wydzia l Fizyki i Informatyki Stosowanej

Praca in ˙zynierska

Daniel Osielczak

kierunek studi´

ow: fizyka techniczna

specjalno´

s´

c: fizyka komputerowa

Generacja liczb przypadkowych przy

u ˙zyciu metody Metropolisa

Opiekun: prof. dr hab. Stanis law Bednarek

Ocena: ................................. Data: ................................. Podpis: .................................

Krak´

ow, stycze´

n 2008

O´swiadczam, ´swiadomy(-a) odpowiedzialno´sci karnej za po´swiadczenie nieprawdy, ˙ze ni-

niejsz

,

a prac

,

e dyplomow

,

a wykona lem(-am) osobi´scie i samodzielnie i nie korzysta lem(-am) ze

´

zr´

ode l innych ni˙z wymienione w pracy.

Skala ocen: (6.0 – celuj

,

aca), 5.0 – bardzo dobra, 4.5 – plus dobra, 4.0 – dobra, 3.5 – plus dostateczna, 3.0 – dostateczna, 2.0 – niedostateczna

Spis tre´

sci

1

Wst

,

ep

4

1.1

Cele pracy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

1.2

Podstawowe poj

,

ecia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

2

Opis metod

5

2.1

Metoda eliminacji von Numanna . . . . . . . . . . . . . . . . . . . . . . . . . . .

5

2.2

Metoda Metropolisa

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

3

Opis test´

ow dla generator´

ow liczb losowych

9

3.1

Testy zgodno´sci rozk ladu . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

9

3.1.1

Test chi-kwadrat . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

9

3.1.2

Test Ko lmogorowa . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

10

3.2

Testy nizale˙zno´sci . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

10

4

Komputerowa realizacja zagadnienia

12

5

Otrzymane wyniki

14

5.1

Jednowymiarowy rozk lad normalny . . . . . . . . . . . . . . . . . . . . . . . . .

14

5.2

Jednowymiarowy rozk lad Erlanga . . . . . . . . . . . . . . . . . . . . . . . . . .

18

5.3

Dwuwymiarowy rozk lad normalny . . . . . . . . . . . . . . . . . . . . . . . . . .

21

6

Podsumowanie

25

Spis rysunk´

ow

26

Spis tabel

27

Lteratura

28

3

1

Wst

,

ep

1.1

Cele pracy

Celem pracy by lo zbadanie w la´sciwo´sci generator´

ow liczb pseudolosowych dzia laj

,

acych w opar-

ciu o algorytmy Metropolisa i von Neumanna. Dla obu metod badania prowadzone by ly na
rozk ladach jedno- i dwuwymiarowych. Wykonane zosta ly tak˙ze por´

ownania efektywno´sci obu

metod we wszystkich wymienionych przypadkach.

1.2

Podstawowe poj

,

ecia

Liczbami losowymi nazywamy ci

,

ag liczb, pomi

,

edzy kt´

orymi nie zachodzi zale˙zno´s´

c funkcyjna,

tzn. pojawienie si

,

e danej liczby w tym ci

,

agu mo˙zna przewidzie´

c tylko z pewnym prawdopo-

dobie´

nstwem. W dalszej cz

,

e´sci pracy b

,

edziemy tak nazywa´

c ci

,

agi liczb uzyskane przy pomocy

programowego generatora liczb losowych (nazywane te˙z liczbami pseudolosowymi gdy˙z faktycz-
nie nie s

,

a dzie lem przypadku, lecz wynikiem procedur matematycznych).[2]

Generatorem liczb losowych (pseudolosowych) b

,

edziemy nazywa´

c program komputerowy

realizuj

,

acy matematyczn

,

a procedur

,

e (algorytm) tworzenia kolejnych liczb poprzez odpowied-

nie dzia lania na liczbie poprzedzaj

,

acej. Dzia lanie generatora mo˙zemy uzna´

c za w la´sciwe je˙zeli

otrzymane przy jego pomocy ci

,

agi liczb przejd

,

a pomy´slnie odpowiednie testy statystyczne.

Testami statystycznymi dla generator´

ow liczb losowych b

,

edziemy nazywa´

c seri

,

e ope-

racji matematycznych wykonanych na ci

,

agach liczb otrzymanych z tych generator´

ow. Celem

ka˙zdego testu losowo´sci jest sprawdzenie, czy badana sekwencja wykazuje cechy prawdziwej
sekwencji losowej lub sprawdzenie, czy nie ma pewnych defekt´

ow statystycznych.

4

2

Opis metod

W rozdziale tym kr´

otko opisane s

,

a metody bed

,

ace obiektem tej pracy. Podane s

,

a ich algorytmy

oraz kilka dodatkowych informacji og´

olnych dotycz

,

acych ka˙zdej z nich.

2.1

Metoda eliminacji von Numanna

Metoda ta pozwala na otrzymanie liczb losowych podlegaj

,

acych penwemu, z g´

ory ustalo-

nemu rozk ladowi. Nale˙zy ona do grupy metod obliczeniowych znanych jako ”Metody Monte
Carlo”, By la to jedna z pierwszych metod pozwalaj

,

acych na generowanie zmiennych o r´

o˙znych

za lo˙zonych rozk ladach. Zosta la ona po raz pierwszy przedstawiona w pracy [3].

Algorytm w najprostszej wersji w przypadku jednowymiarowym przedstawia si

,

e nast

,

epuj

,

aco[1]:

Zmienna losowa X ma rozk lad o g

,

esto´sci r´

ownej f (x) na przedziale (a, b), gdzie −∞ < a < b <

+∞ oraz r´

ownej zeru poza tym przedzia lem. Ponadto f (x) ≤ c dla pewnej sta lej c < +∞.

1. Generuje si

,

e punkt (R

1

, R

2

), przy czym R

1

jest zmienn

,

a losow

,

a o rozk ladzie r´

ownomiernym

na przedziale (a, b), R

2

jest zmienn

,

a losow

,

a o rozk ladzie r´

ownomiernym na przedziale

(0, c), oraz zmienne losowe R

1

i R

2

s

,

a niezale˙zne. Korzystamy przy tym z faktu, ˙ze

je˙zeli zmienna losowa R ma rozk lad r´

ownomierny na przedziale (0, 1), to zmienna losowa

(b − a)R + a ma rozk lad r´

ownomierny na przedziale (a, b).

2. Je˙zeli R

2

≤ f (R

1

), przyjmuje si

,

e X = R

1

, w przypadku przeciwnym par

,

e (R

1

, R

2

)

,

eliminuje si

,

e”i powtarza obliczenia od wygenerowania nowej pary zgodnie z punktem 1.

Otrzymana w ten spos´

ob zmienna losowa X ma rozk lad o g

,

esto´sci f .

Mo˙zna to udowodni´

c w nast

,

epuj

,

acy spos´

ob[1]:

latwo zauwa˙zy´

c, ˙ze zmienna losowa (R

1

, R

2

) ma rozk lad r´

ownomierny w prostok

,

acie o podstawie

Ω i wysoko´sci c, gdzie Ω jest odcinkiem (a, b) (w przypaku zmiennych losowych m-wymiarowych
jest on obszarem m-wymiarowym, w kt´

orem mo˙zna okre´sli´

c rozk lad r´

ownomierny). g

,

eso´s´

c tego

rozk ladu jest wi

,

ec r´

owna 1/c|Ω|, gdzie |Ω| jest d lugo´sci

,

a odcinka Ω.

Mamy wi

,

ec:

P{X ≤ x} = P{R

1

≤ x|R

2

≤ f (R

1

)} =

P{R

1

≤ x, R

2

≤ f (R

1

)}

P{R

2

≤ f (R

1

)}

Ale:

P{R

2

≤ f (R

1

)} =

Z

Ω

 Z

f (R

1

)

0

dR

2

c



dR

1

|Ω|

=

1

c|Ω|

Z

Ω

f (R

1

)dR

1

=

1

c|Ω|

oraz

P{R

1

≤ x, R

2

≤ f (R

1

)} =

1

c|Ω|

Z

x

0

dR

1

Z

f (R

1

)

0

dR

2

=

1

c|Ω|

Z

x

0

f (R

1

)dR

1

Otrzymali´smy wi

,

ec:

P{X ≤ x} =

Z

x

0

f (R

1

)dR

1

5

Rysunek 1: Metoda von Neumanna

Sk

,

ad bezpo´srednio wynika, ˙ze g

,

esto´sci

,

a zmiennej losowej X jest f (x).

Dow´

od ten mo˙zna bardzo latwo rozsze˙zy´

c na zmienne losowe m-wymiarowe otrzymywane t

,

a

metod

,

a.

Algorytm elimiancji von Neumanna dla m-wymiarowych zmiennych losowych jest analogiczny
do podanego wcze´sniej jednowymiarowego. W przypadku rozk lad´

ow m-wymiarowych pewien

dodatkowy problem mo˙ze stanowi´

c losowanie wed lug rozk ladu r´

ownomiernego na zbiorze Ω.

Trudno jest sformu lowa´

c jakie´s og´

olne regu ly przeprowadzenia takiego losowania. Gdy Ω jest

ograniczonym obszarem w m-wymiarowej przestrzeni liczb rzeczywistych, najprostsze podej´scie
polega na zbudowaniu takiego m-wymiarowego przedzia lu
{(x

1

, x

2

, . . . , x

m

) : a

1

≤ x

1

≤ b

1

, a

2

≤ x

2

≤ b

2

, . . . , a

m

≤ x

m

≤ b

m

, }

wewn

,

atrz kt´

orego znajduje si

,

e Ω. Wtedy wsp´

o lrz

,

edne X

1

, X

2

, . . . , X

m

losuje si

,

e niezale˙znie

wed lug rozk ladu r´

ownomiernego na przedziale, odpowiednio (a

1

, b

1

), (a

2

, b

2

), . . . , (a

m

, b

m

) i punkt

(X

1

, X

2

, . . . , X

m

) akceptuje si

,

e do oblicze´

n gdy nale˙zy do Ω, a w przypadku przeciwnym od-

rzuca i powtarza losowanie.

Ograniczenie przedzia lu od g´

ory sta l

,

a powoduje jednak s lab

,

a jako´s´

c generowanych rozk lad´

ow.

Dlatego cz

,

e´sciej wykorzystywan

,

a jest metoda [6]: Szukamy zmiennej losowej na przdziale (a, b)

podlegaj

,

acej rozk ladowi f (x).

Algorytm omawianej metody jest w´

owczas nast

,

epuj

,

acy:

1. W przedziale (a, b) okre´slamy taki rozk lad g(x), ˙ze dla ka˙zdego x ∈ (a, b) : g(x) > f (x).

2. Stosuj

,

ac dowolny generator o rozk ladzie r´

ownomiernym w przedziale (a, b) losujemy nowy

punkt x

n

, a nast

,

epnie wyznaczamy w

n

= f (xn)/g(xn).

3. Stosuj

,

ac dowolny generator o rozk ladzie r´

ownomiernym w przedziale (0, 1) losujemy liczb

,

e

6

losow

,

a r

n

.

4. Je˙zeli r

n

< w

n

, punkt x

n

przyjmujemy, je˙zeli nie - odrzucamy.

Problemem mo˙ze tu by´

c w niekt´

orych przypadkach wyznaczenie g(x), gdy˙z nie mo˙ze ona by´

c

znacznie wi

,

eksza od f (x) na tym przedziale ale jednocze´snie nie mo˙ze by´

c mniejsza.

Uog´

olnieniem tej metody pozwalaj

,

acym na konstruowanie prostych algorytm´

ow generowania

r´

o˙znych zmiennych losowych jest post

,

epowanie:

Niech zmienna losowa X ma rozk lad o g

,

esto´sci p(x) = a

1

f

1

(x)g

1

(x)

gdzie f

1

jest pewn

,

a g

,

esto´sci

,

a prawdopodobie´

nstwa, funkcja g

1

spe lnia warunek 0 ≤ g

1

(x) ≤ 1

dla ka˙zdego x oraz a

1

jest sta l

,

a. Otrzymujemy nast

,

epuj

,

acy algorytm:

1. Wygenerowa´

c zmienn

,

a y wed lug rozk ladu o g

,

esto´sci f

1

(y)

2. Obliczy´

c g

1

(y)

3. Wygenerowa´

c R wed lug rozk ladu r´

ownomiernego na przedziale (0, 1)

4. Je˙zeli R ≤ g

1

(y), przyj

,

a´

c X = y, w przypadku przeciwnym odrzuci´

c wygenerowane liczby

i powt´

orzy´

c obliczenia od pkt 1.

Poniewa˙z dla tak skonstruowanej zmiennej losowej X mamy

P{X ≤ x} = P{y ≤ x|R ≤ g

1

(y)} =

Z

R≤g

1

(y)

Z

y≤x

f

1

(y)dydR

Z

R≤g

1

(y)

Z

f

1

(y)dydR

= a

1

Z

x

−∞

 Z

g

1

(x)

0

dR



f

1

(y)dy =

a

1

Z

x

−∞

f

1

(y)g

1

(y)dy =

Z

x

−∞

p(y)dy

czyli prowadzi to do wygenerowania zmiennej losowej o rozk ladzie z g

,

esto´sci

,

a p(x).

Rozsze˙zenie powy˙zszej metody na przypadek wielowymiarowy nie stwarza zbyt du˙zo pro-
blem´

ow.

Warto zaznaczy´

c, ˙ze algorytm ten nie jest sam w sobie generatorem, a jedynie filtrem ko-

rzystaj

,

acym z liczb losowych wygenerowanych przy u˙zyciu innych generator´

ow o rozk ladzie

r´

ownomiernym.

Tak wi

,

ec du˙zy wp lyw na jako´s´

c (zw laszcza losowo´s´

c”) otrzymanych liczb

b

,

edzie mia l sam generator o rozk ladzie r´

ownomiernym. Dla s labych generator´

ow odrzucone

zostanie zbyt du˙zo liczb pr´

obnych i wydajno´s´

c metody zmniejszy si

,

e znacznie.

Kolejnym problemem tej metody jest radzenie sobie z rozk ladami (zw laszcza wielowymiaro-
wymi), kt´

ore posiadaj

,

a silne lokalne maksima. Rozk lad w takich wypadkach odbiega zazwyczaj

znacznie od za lo˙zonego.

2.2

Metoda Metropolisa

Metoda ta, podobnie jak metoda eliminacji von Neumanna, pozwala na otrzymanie liczb lo-
sowych podlegaj

,

acych pewnemu, z g´

ory ustalonemu rozk ladowi. Zosta la przedstawiona po raz

pierwszy w pracy [4].

Jest to najpopularniejszy obecnie algorytm z grupy ”Metod Monte Carlo”. Ide

,

a jest b l

,

adzenie

7

przypadkowe, kt´

ore przez odrzucanie cz

,

e´sci z wylosowanych liczb prowadzi do otrzymania

zmiennych wed lug zadanego rozk ladu.

Jest wi

,

ec, tak jak poprzednia metoda, filtrem ko-

rzystaj

,

acym z liczb losowych wygenerowanych przy u˙zyciu innych generator´

ow o rozk ladzie

r´

ownomiernym.

Metoda ta jest powszechnie u˙zywana w przypadku rozk lad´

ow, dla kt´

orych metoda von Neu-

manna okazuje si

,

e nieefektywna.

Algorytm przedstawiony zosta l dla problemu jednowymiarowego ale mo˙zna go latwo rozsze˙zy´

c

na m wymiar´

ow post

,

epuj

,

ac podobnie jak w przypadku metody von Neumanna.

x

1

, x

2

, . . . , x

n

jest zbiorem liczb losowych o zadanym rozk ladzie f (x). Kolejne jego elementy

otrzymujemy w nast

,

epuj

,

acy spos´

ob [6]:

1. Wybieramy punkt startowy x

i

. Mo˙zemy to zrobi´

c losowo, ale zazwyczaj wybiera si

,

e go

poprzez wykonanie kilku pocz

,

atkowych krok´

ow algorytmu i wybranie po nich tego punktu

dla kt´

orego warto´s´

c f (x

i

) jest dostatecznie du˙za [5]

2. Wyznaczamy nowy punkt x

t

poprzez wylosowanie go z otoczenia x

i

: x

i

− ∆x ≤ x

t

≤

x

i

+∆x, oraz warto´s´

c f (x

t

) w tym punkcie. ∆x dobieramy w ten spos´

ob aby oko lo po lowy

pr´

obnych liczb zosta lo zaakceptowanych.

3. Je˙zeli warto´s´

c f (x

t

) ≥ f (x

i

) to akceptujemy nowy punkt i przyjmujemy x

i+1

= x

t

. W

wypadku f (x

t

) < f (x

i

) generujemy wed lug rozk ladu jednorodnego liczb

,

e losow

,

a r z

przedzia lu (0, 1):

• je˙zeli r < f (x

t

)/f (x

i

), akceptujemy nowy punkt i przyjmujemy x

i+1

= x

t

• je˙zeli r ≥ f (x

t

)/f (x

i

),odrzucamy nowy punkt i przyjmujemy x

i+1

= x

i

4. x

i+1

staje si

,

e nowym punktem startowym (x

i

= x

i+1

)

Po wykonaniu odpowiednio du˙zej liczby krok´

ow procedura doprowadza do zbioru punkt´

ow x

zgodnego z za lo˙zonym rozk ladem f (x).

Algorytm Metropolis definiuje zbie˙zny la´

ncuch Markowa, kt´

orego granica jest szukanym rozk ladem

prawdopodobie´

nstwa. Jednak nie we wszystkich problemach zbie˙zno´s´

c ta jest wystarczaj

,

aco

szybka. W niekt´

orych wypadkach liczby krok´

ow algorytmu mo˙ze by´

c ograniczona wielomia-

nowo, w innych mo˙zna dowie´s´

c, ˙ze takie ograniczenie nie istnieje.

8

3

Opis test´

ow dla generator´

ow liczb losowych

Testowanie generator´

ow liczb (pseudo)losowych polega na badaniu ci

,

ag´

ow liczb produkowa-

nych przez ten generator za pomoc

,

a odpowiednich test´

ow statystycznych. Kolejn

,

a liczb

,

e X

n

produkowan

,

a przez ten generator traktujemy jako zmienn

,

a losow

,

a. Weryfikacja, czy generator

produkuje liczby losowe o ˙z

,

adanym rozk ladzie prawdopodobie´

nstwa, sprowadza si

,

e wi

,

ec do we-

ryfikacji, czy ci

,

ag X

0

, X

1

, . . . , X

N −1

mo˙ze by´

c traktowany jako N -elementowa pr´

obka prosta z

okre´slonej populacji, czyli, czy jest on ci

,

agiem niezale˙znych zmiennych losowych o jednakowym

rozk ladzie prawdopodobie´

nstwa. Zaakceptowanie generatora na podstawie takiej analizy jest

jednak niemo˙zliwe. Je˙zeli jednak generator pomy´slnie przeszed l wystarczaj

,

ac

,

a liczb

,

e test´

ow na

poziomie istotno´sci α, to mo˙zemy go zaakceptowa´

c, opieraj

,

ac si

,

e na zaufaniu, i˙z przejdzie on

pomy´slnie tak˙ze inne testy.[1]

Testy statystyczne przeznaczone do weryfikacji generator´

ow liczb losowych mo˙zna podzieli´

c na

dwie grupy:

1. testy zgodno´sci rozk ladu

2. testy niezale˙zno´sci (losowo´sci pr´

obki)

3.1

Testy zgodno´

sci rozk ladu

Za pomoc

,

a tych test´

ow weryfikujemy hipotez

,

e, ˙ze pr´

obka X

0

, X

1

, . . . , X

N −1

pochodzi z populacji

o okre´slonym rozk ladzie prawdopodobie´

nstwa. W lasno´sci takiej pr´

obki to np[1]:

1. Je˙zeli X

0

, X

1

, . . . , X

N −1

jest pr´

obk

,

a prost

,

a populacji o rozk ladzie r´

ownomiernym na prze-

dziale (0, 1), to warto´s´

c oczekiwana ´sredniej arytmetycznej z pr´

obki X

N

=

P

X

i

/N = 1/2.

Oznacza to, ˙ze ´srednia liczb wyprodukowanych przez generator nie powinna si

,

e znacz

,

aco

r´

o˙zni´

c od 1/2.

2. Je˙zeli X

0

, X

1

, . . . , X

N −1

jest pr´

obk

,

a prost

,

a populacji o rozk ladzie r´

ownomiernym na prze-

dziale (0, 1), to warto´s´

c oczekiwana ´sredniej kwadratu X

2

N

=

P

X

2

i

/N = 1/3. Oznacza

to, ˙ze ´srednia liczb wyprodukowanych przez generator nie powinna si

,

e znacz

,

aco r´

o˙zni´

c od

1/3.

3. Podzielmy przedzia l (0, 1) punktami 0 = a

0

< a

1

< . . . < a

s

= 1 na podprzedzia ly.

Je˙zeli generator produkuje liczby lkosowe orozk ladzie r´

ownomiernym na przedziale (0, 1),

to frakcja tych liczb w ci

,

agu x

0

, x

1

, . . . , x

N −1

, kt´

ore nale˙z

,

a do przedzia lu (a

j

, a

j+1

), nie

powinna r´

o˙zni´

c si

,

e istotnie od liczby a

j+1

− a

j

.

Najpopularniejszymi testami tego typu s

,

a mi

,

edzy innymi[7]:

3.1.1

Test chi-kwadrat

Jest to jeden z najbardziej popularnych test´

ow zgodno´sci dla generator´

ow liczb losowych. Jego

stablicowane warto´sci krytyczne mo˙zna znale´

z´

c w tablicach matematycznych, podr

,

ecznikach

czy matematycznych programach komputerowych.
Hipoteza, kt´

or

,

a sprawdzamy: zmienna losowa X ma rozk lad prawdopodobie´

nstwa o dystrybu-

ancie F .
Uzyskane wyniki grupujemy w k klas. Definiujemy nast

,

epuj

,

ac

,

a statystyk

,

e:

9

χ

2
k−1

=

k

X

i=1

(n

i

− np

i

)

2

np

i

gdzie n

i

to liczba takich element´

ow X ci

,

agu X

1

, X

2

, . . . , X

n

spe lniaj

,

acych warunek a

i−1

< X ≤

a

i

(a

i

to miejsca podzia lu zbioru warto´sci zmiennej losowej X na poszczeg´

olne klasy). p

i

to

oczekiwane prawdopodobie´

nstwo danej klasy.

Je˙zeli dobierzemy rozbicie tak, by p

i

= 1/k , w´

owczas statystyka ta wygl

,

ada w nast

,

epuj

,

acy

spos´

ob: χ

2
k−1

=

k

n

k

X

i=1

n

2
i

− n

Wygodnie jest przyj

,

a´

c k r´

owne oko lo 10, w´

owczas latwo obliczy´

c poziomy krytyczne F (χ

2
k−1

)

testu.

3.1.2

Test Ko lmogorowa

Weryfikowana hipoteza: zmienna losowa X ma rozk lad o danej ci

,

ag lej dystrybuancie F , a

statystyka testu jest oparta na r´

o˙znicy pomi

,

edzy hipotetyczn

,

a dystrybuant

,

a F a dystrybuant

,

a

rzeczywist

,

a F

n

, opart

,

a na pr´

obie X

1

, X

2

, . . . , X

n

.

Statystyk

,

a testow

,

a jest:

D

n

=

sup

− inf<x<inf

|F

n

(x) − F (x)|

Dystrybuanta rzeczywista jest okre´slona wzorem:

F

n

(x) = n

−1

n

X

j=1

1

(− inf,x]

(X

j

)

Je˙zeli pr´

oba X

1

, X

2

, . . . , X

n

powsta la z rozk ladu okre´slonego przez F , w´

owczas warto´s´

c D

n

d

,

a˙zy do zera z prawdopodobie´

nstwem r´

ownym 1.

Je˙zeli ta statystyka osi

,

aga du˙ze warto´sci, nale˙zy zastanowi´

c si

,

e nad odrzuceniem weryfikowanej

hipotezy.
Aby obliczy´

c statystyk

,

e D

n

z pr´

oby X

1

, X

2

, . . . , X

n

nale˙zy zauwa˙zy´

c, ˙ze D

n

= sup

− inf<x<inf

|F

n

(x)−

F (x)| jest osi

,

agane w pewnym punkcie, w kt´

orym nast

,

epuje skok dystrybuanty F

n

. Poniewa˙z

D

n

jest sta le przy monotonicznych przekszta lceniach x, w´

owczas:

D

+

n

= max

1≤i≤n

(

i

n

− F (X

i:n

))

D

−

n

= max

1≤i≤n

(F (X

i:n

) −

i − 1

n

)

D

n

= max{D

+

n

, D

−

n

}

zapis X

i:n

oznacza i-t

,

a statystyk

,

e pozycyjn

,

a z danej pr´

oby.

Warto´sci krytyczne dla tego testu r´

ownie˙z s

,

a stablicowane.

3.2

Testy nizale ˙zno´

sci

S lu˙z

,

a do weryfikacji tezy, ˙ze ci

,

ag X

0

, X

1

, . . . , X

N −1

jest pr´

obk

,

a prost

,

a, czyli, czy kolejno pro-

dukowane przez generator liczby mog

,

a by´

c traktowane jako realizacje niezale˙znych zmiennych

losowych. Wyr´

o˙zniamy trzy grupy test´

ow[1]:

1. testy zgodno´sci rozk lad´

ow wielowymiarowych

10

2. testy autokorelacji

3. testy uporz

,

adkowania pr´

obki (kombinatoryczne)

Testy zgodno´

sci rozk lad´

ow wielowymiarowych opieraj

,

a si

,

e zazwyczaj na spostrze˙zeniu,

˙ze je˙zeli X

0

, X

1

, . . . , jest ci

,

agiem liczb z generatora liczb losowych o rozk ladzie r´

ownomiernym

na przedziale (0, 1), to k-wymiarowe zmienne losowe

(X

jk

, X

jk+1

, . . . , X

(j+1)(k−1)

)j = 0, 1, 2, . . .

s

,

a niezale˙zne i maj

,

a jednakowy rozk lad r´

ownomierny na kostce jednostkowej (0, 1)

k

:

(0, 1)

k

= {x = (x

1

, x

2

, . . . , x

k

) : 0 < x

j

< 1, j = 1, 2, . . . , k}

Testy autokorelacji opieraj

,

a si

,

e na spostrze˙zeniu, ˙ze je˙zeli X

0

, X

1

, . . . , jest ci

,

agiem liczb z

generatora liczb losowych o rozk ladzie r´

ownomiernym na przedziale (0, 1), to wsp´

o lczynnik

korelacji %(X

n

, X

n+r

) nie zale˙zy od n i jest r´

owny zeru dla wszystkich r ≥ 1 (ci

,

ag o takiej

w lasno´sci nazywa si

,

e zazwyczaj ci

,

agiem bia lym). Testowanie hipotezy, ˙ze tak rzeczywi´scie jest,

polega na sprawdzeniu czy wsp´

olczynniki korelacji

ˆ

%

r

=

1

(N − r)

N −r−1

X

n=0



x

n

−

1

2



x

n+r

−

1

2



obliczone dla zaobserwanego ci

,

agu x

0

, x

1

, . . . , x

N −1

liczb rozpatrywanego generatora r´

o˙zni

,

a si

,

e

istotnie od zera.[1]

Testy uporz

,

adkowania pr´

obki to wszelkiego rodzaju testy serii, testy statystyk pozycyjnych,

testy kombinatoryczne, takie jak np. test pokerowy itd. Testu te maj

,

a za zadanie wykrycie

ewentualnych tendencji w ci

,

agach liczb produkowanych przez generator, co prowadzi loby do

jego dyskwalifikacji. Nie b

,

edziemy si

,

e tutaj nimi dok ladnie zajmowa´

c. wi

,

ecej informacji mo˙zna

znale´

z´

c w pozycjach [1], [7] i [8]

11

4

Komputerowa realizacja zagadnienia

Zagadnienie zosta lo zaprogramowane w j

,

ezyku C++[10], wykorzystane zosta ly standardowe bi-

blioteki iostream s lu˙z

,

aca do obs lugi standardowego wej´scia/wyj´scia, fstream s lu˙z

,

aca do bs lugi

zapisu i odczytu plik´

ow, time.h potrzebna do ustanowienia losowego seed dla generatora biblio-

tecznego liczb losowych, oraz math.h zawieraj

,

aca definicj

,

e procedur matematycznych takich

jak sqrt() czy exp().
Pierwszym problemem napotkanym podczas realizacji zagadnienie by l dob´

or generatora liczb

losowych. Potrzebny by l generator daj

,

acy liczby o rozk ladzie r´

ownomiernym na przedziale

[0, 1). Dodatkowo powinien on by´

c na tyle szybki aby nie obci

,

a˙za´

c nadmiernie procedury sa-

mym losowaniem liczb. Najprostszym wyborem by la oczywi´scie funkcja biblioteczna j

,

ezyka

C++ rand(). Pomimo krytyki tej funkcji zamieszczonej w [12] okaza la si

,

e ona szybsz

,

a i daj

,

ac

,

a

lepsze rozk lady jednorodne ni˙z testowany generator mieszany postaci [14]:

ξ

i+1

=

(aξ

i

+ b)mod m

m − 1

zw laszcza dla ci

,

ag´

ow liczb rz

,

edu 1000000 ( z powodu kr´

otkiego okresu ).

Najlepszym rozwi

,

azaniem wydaje si

,

e jednak generator oparty na algorymie typu Mersenne

Twister, kt´

orego dok ladny opis znale´

z´

c mo˙zna np. w [13]. Pewnym problemem w tym wypadku

by la implementacja dlatego zdecydowano si

,

e na gotow

,

a bibliotek

,

e w j

,

ezyku C++ dost

,

epn

,

a pod

adresem [15]. Jest ona dost

,

epna na licencji GNU.

Wykresy 2 (generator mieszany), 3 (funkcja rand()), 4 ( Marsenne Twister) pokazuj

,

a rozk lad

wylosowanych liczb z przedzia lu [0, 1) dla ka˙zdego z generator´

ow. W tabeli zebrano ´srednie z

1000000 liczb dla kolejnych generator´

ow w tym przedziale (aby rozk lad by l jednorodny powinny

by´

c bliskie 1/2).

Tablica 1: ´

Srednia liczb na przedziale [0,1)

Typ generatora

´

srednia

mieszany

0.500001

rand()

0.499255

Marsenne Twister

0.499387

Jak wida´

c wszystkie z testowanych generator´

ow maj

,

a rozk lad wystarczaj

,

aco bliski 1/2. Na

Marsenne Twister zdecydowano si

,

e g l´

ownie z powodu jego d lugiego okresu.

Drug

,

a kwesti

,

a by l podzia l zakresu otrzymanych liczb na fragmenty. By lo to konieczne aby

otrzyma´

c rozk lad g

,

esto´sci prawdopodobie´

nstwa. Zdecydowano si

,

e na dzielenie w stosunku 1 : 10

tj. ka˙zdy odcinek o d lugo´sci r´

ownej jedno´sci podzielony zosta l na 10 pododcink´

ow. Dok ladno´s´

c

takiego podzia lu jest wystarczaj

,

aca dla test´

ow statystycznych a jednocze´snie niezbyt obci

,

a˙za

obliczenia.
Kolejn

,

a kwesti

,

a by la wizualizacja wynik´

ow kt´

orej dokonano za pomoc

,

a programu gnuplot

wchodz

,

acego w sk lad standardowych narz

,

edzi w systemie Linux. Daje on du˙ze mo˙zliwo´sci

zar´

owno je´sli chodzi o wykresy dwuwymiarowe, jak i wykresy 3D.

12

Rysunek 2: Generator mieszany

Rysunek 3: Funkcja rand()

Rysunek 4: Marsenne Twister

13

5

Otrzymane wyniki

5.1

Jednowymiarowy rozk lad normalny

Pierwszym zrealizowanym zadaniem by la generacja standardowego rozk ladu normalnego ( ozna-
czanego czasami N (0, 1) ) tj. rozk ladu normalnedo z warto´sci

,

a oczekiwan

,

a µ = 0 oraz wa-

riancj

,

a σ

2

= 1. Rozk lad ten opisuje funkcja g

,

esto´sci prawdopodobie´

nstwa postaci f (x) =

1

√

2π

exp (

−x

2

2

) .

Generacja wykonana zosta la przy u˙zyciu obu metod ( von Neumanna i Metropolisa). Dla
ka˙zdego przypadku generowane by lo 1000000 liczb. Jest to ilo´s´

c optymalna, poniewa˙z oblicze-

nia wykonywane s

,

a dla tej niej szybko ( czas ok 10 sekund ), a jednocze´snie ilo´s´

c liczb jest

wystarczaj

,

aca do przeprowadzenia test´

ow statystycznych.

Wykresy 5 i 6 prezentuj

,

a otrzymane liczby.

Dla obu metod sprawdzona zosta la ilo´s´

c liczb odrzuconych spo´sr´

od wylosowanych.

Dane dla metody Metropolisa zosta ly uzyskane dla warto´sci ∆ = 1.9. Jest to warto´s´

c najbli˙zsza

optymalnej ( powinno si

,

e j

,

a dobiera´

c w ten spos´

ob aby ok. 2/3 warto´sci by lo akceptowanych).

Natomiast dla metody von Neumanna posta´

c funkcji g(x) zosta la wybrana jako g(x) = 0.4, gdy˙z

jest to g´

orna granica warto´sci osi

,

aganych przez f (x) =

1

√

2π

exp (

−x

2

2

) . Ograniczenie takie oka-

zuje si

,

e dawa´

c wystarczaj

,

aco dobre wyniki, w przeciwie´

nstwie do funkcji g(x) =

s

2e

π

1

2

exp−|x|

sugerowanej w pracy [9]. Funkcja ta dawa la niezadowalaj

,

ace wyniki, gdy˙z prawdopodobie´

nstwo

osi

,

aga lo tam maksima dla warto´sci |x| = 1 co wida´

c na wykresie 7.

Dla obu metod wykonane zosta lo po 4 losowania, a nast

,

epnie policzona zosta la ´srednia. Ilo´sci

odrzuconych liczb dla ka˙zdego losowania oraz ´srednie podane zosta ly w tabeli 2.

Tablica 2: Ilo´s´

c odrzuconych liczb na 1000000 wylosowanych dla N(0,1) w 2D

Metoda

losowanie 1

losowanie 2

losowanie 3

losowanie 4

´

srednia

von Neumanna

2197421

2199050

2200888

2197666

2198756.25

Metropolisa

352338

353227

352620

352806

352747.75

Nast

,

epnie sprawdzono efektywno´s´

c metod, rozumian

,

a jako stosunek liczby zaakceptowanych

zmiennych losowych do liczby wszystkich wygenerowanych zmiennych losowych (tj. zaakcep-
towanych + odrzuconych). Otrzymane wyniki dla 4 losowa´

n, oraz ich ´srednia zebrane zosta ly

w tabeli 3. Wyra´

znie widoczna jest tu przewaga metody Metropolisa, kt´

ora wymaga wygene-

Tablica 3: Stosunek zaakceptowane/wylosowane dla N(0,1) w 2D

Metoda

losowanie 1

losowanie 2

losowanie 3

losowanie 4

´

srednia

von Neumanna

31.2752059%

31.2592801%

31.2413305%

31.2728096%

31.262156525%

Metropolisa

73.9460105%

73.8974318%

73.930594%

73.9204291%

73.92361635%

rowania ok 2,5 razy mniejszej ilo´sci liczb. Przek lada si

,

e to na szybko´s´

c dzia lania algorytmu.

R´

o˙znica jest zauwa˙zalna, ale znacznie mniejsza ni˙z ten stosunek, mianowicie ok 20% na korzy´s´

c

14

Rysunek 5: Liczby otrzymane metod

,

a von Neumanna dla rozk ladu N(0,1) w 2D

Rysunek 6: Liczby otrzymane metod

,

a Metropolisa dla rozk ladu N(0,1) w 2D

15

Rysunek 7: Rozk lad liczb otrzymanych metod

,

a von Neumanna dla g(x) = 1/2 * sqrt(2e/π) *

exp(-|x|)

metody Metropolisa ( poniewa˙z metoda von Neumanna wymaga mniejszej ilo´sci por´

owna´

n).

Nie zmienia si

,

e ona znacznie przy zmianie ilo´sci generowanych liczb.

Nast

,

epnie przyporz

,

adkowano liczby do zakres´

ow o szeroko´sci 1/10. Dzi

,

eki temu otrzymano

rozk lady g

,

esto´sci pokazane na wykresach 8 i 9. Na wykresach umieszczono r´

ownie˙z funkcj

,

e

g

,

esto´sci prawdopodobie´

nstwa dla N(0,1), czyli f (x) =

1

√

2π

exp (

−x

2

2

) pomno˙zon

,

a 100000 razy.

Jest to liczba 10-krotnie mniejsza od ilo´sci liczb wylosowanych, a jednak zgadza si

,

e ona z

g

,

esto´sci

,

a. Dzieje si

,

e tak poniewa˙z f (x) jest funkcj

,

a ci

,

ag l

,

a i aby por´

ownywa´

c j

,

a z dyskretnym

rozk ladem nale˙zy przemno˙zy´

c j

,

a przez podzia l zastosowany do owego rozk lady wynosz

,

acy tu

1 : 10.
Ostatnim zadaniem wykonanym dla tego rozk ladu by lo przeprowadzenie test´

ow statystycznych

otrzymanych liczb pod k

,

atem zgodno´sci z rozk ladem N(0,1). Dokonano dw´

och niezale˙znych

test´

ow, mianowicie testu chi-kwadrat opisanego we wst

,

epie oraz testu Ko lmogorowa-Smirnowa,

bed

,

acego jedn

,

a z wersji opisanego testu Ko lmogorowa. Warto´sci dystrybunaty rozk ladu nor-

malnego oraz warto´sci krytyczne dla test´

ow zaczerpni

,

ete zosta ly ze strony [16].

Otrzymane wyniki dla testu chi-kwadrat prezentuje tabela 4. Otrzymane warto´sci por´

ownano

Tablica 4: Wyniki testu chi-kwadrat dla obu metod

Metoda

losowanie 1

losowanie 2

losowanie 3

losowanie 4

´

srednia

von Neumanna

89.6059

79.9145

85.9973

87.4262

85.735975

Metropolisa

85.883

73.597

84.162

94.5955

84.559375

z warto´sci

,

a krytyczn

,

a tego testu dla poziomu istotno´sci α = 0.05 oraz ilo´sci stopni swobody

df = 79. Wynosi la ona χ

2
α,df

= 100.7486, co jest warto´sci

,

a wi

,

eksz

,

a od otrzymanych, wi

,

ec test

ten nie daje nam podstaw do odrzucenia hipotezy, i˙z otrzymane liczby maj

,

a rozk lad zgodny ze

standardowym rozk ladem normalnym N(0,1). Dla otrzymanych warto´sci hipotez

,

e t

,

a da loby si

,

e

odrzuci´

c dla warto´sci α ≥ 0.283, co jest ju˙z warto´sci

,

a bardzo wysok

,

a dla tak du˙zej ilo´sci liczb.

16

Rysunek 8: Rozk lad liczb otrzymanych metod

,

a von Neumanna dla rozk ladu N(0,1) w 2D

Rysunek 9: Rozk lad liczb otrzymanych metod

,

a Metropolisa dla rozk ladu N(0,1) w 2D

17

Drugim z zastosowanych test´

ow by l test Ko lmogorowa-Smirnowa kt´

orego wyniki prezentuje ta-

bela 5. Otrzymane warto´sci por´

ownano z warto´sci

,

a krytyczn

,

a tego testu dla poziomu istotno´sci

Tablica 5: Wyniki testu Ko lmogorowa-Smirnowa dla obu metod

Metoda

losowanie 1

losowanie 2

losowanie 3

losowanie 4

´

srednia

von Neumanna

0.17156

0.17156

0.175934

0.183349

0.17560075

Metropolisa

0.180415

0.178993

0.179055

0.181819

0.1800705

α = 0.05 oraz znanej liczby pomiar´

ow, kt´

ora wynosi la w tym wypadku λ

k

= 0, 152052622. Jest

to warto´s´

c mniejsza od otrzymanych nale˙zy wi

,

ec zmniejszy´

c poziom istotno´sci do α = 0.01, dla

kt´

orego warto´s´

c λ

k

= 0, 18223954. Przy tym poziomie nie mamy ju˙z powod´

ow do odrzucenia

hipotezy, i˙z otrzymane liczby maj

,

a rozk lad zgodny ze standardowym rozk ladem normalnym

N(0,1).

5.2

Jednowymiarowy rozk lad Erlanga

Rozk lad ten zosta l opisany przez A. K. Erlanga do szacowania liczby rozm´

ow telefonicznych,

l

,

aczonych jednocze´snie przez operatora w r

,

ecznej centrali telefonicznej. Obecnie rozk lad ten

znalaz l te˙z zastosowanie w teorii proces´

ow stochastycznych.

Funkcj g

,

esto´sci prawdopodo-

bie´

nstwa ma posta´

c

f(x)=











λ

k

x

k−1

e

−λx

(k − 1)!

dla x ≥ 0

0

dla x < 0

gdzie k > 0 jest pewn

,

a liczb

,

a ca lkowit

,

a zwan

,

a parametrem kszta ltu, a λ > 0 to liczba rzeczy-

wista zwana cz

,

esto´

sci

,

a (odwrotno´s´

c parametru kszta ltu) [2].

Warto´s´

c oczekiwana tego rozk ladu wynosi k/λ.

Rozk lad ten jest specjalnym przypadkiem rozk ladu Gamma. Stosowany tu parametr k jest
liczb

,

a ca lkowit

,

a, w przeciwie´

nstwie do rozk ladu Gamma gdzie ma on posta´

c liczby rzeczywi-

stej. Pozwala to na stosowanie funcji silnia, zamiast funkcji Γ Eulera. Fakt ten znacznie u latwia
implementacj

,

e tego rozk ladu, gdy˙z uwalnia wz´

or od trudnego i generuj

,

acego b l

,

edy ca lkowania

numerycznego.
Tak jak w przypadku poprzedniego podpunktu wygenerowane zosta lo 1000000 liczb. Szybko´sci
dzia lania obu metod nie zmieni ly si

,

e zauwa˙zalnie w stosunku do do wcze´sniejszego zadania.

Pewn

,

a zmian

,

a by la jedynie zmiana parametru ∆, kt´

ory w tym zadaniu ustawiony zosta l na

0.5.
Dla rozk ladu Erlanga przyj

,

eto parametry k = 5 oraz λ = 1.0, poniewa˙z przy tak dobranych

warto´sciach wykres g

,

esto´sci prawdopodobie´

nstwa zanika l (czyli zmierza l do 0) przy warto´sci ok

20. Otrzymane liczby prezentuj

,

a wykresy 10 i 11.

Rozk lad liczb otrzymano poprzez podzielenie przedzia lu [0, 20] na 100 podprzedzia l´

ow co da lo

szeroko´s´

c pojedynczego przedzia lu 0.2. Dlatego te˙z wykres warto´sci funkcji g

,

esto´sci rozk ladu

pomno˙zony jest przez ilo´s´

c wylosowanych liczb a nast

,

epnie przez stosunek podzia lu czyli 1:5.

Rozk lady dla obu metod prezentuj

,

a wykresy 12 i 13.

Nast

,

epnie, tak jak w poprzednim podpunkcie, dla obu metod wykonane zosta lo po 4 losowania

oraz policzona zosta la ´srednia ilo´sci odrzuconych liczb. Wyniki dla ka˙zdego losowania oraz

´srednie podane zosta ly w tabeli 6.

18

Rysunek 10: Liczby otrzymane metod

,

a von Neumanna dla rozk ladu Erlanga

Rysunek 11: Liczby otrzymane metod

,

a Metropolisa dla rozk ladu Erlanga

19

Rysunek 12: Rozk lad liczb otrzymanych metod

,

a von Neumanna dla rozk ladu Erlanga

Rysunek 13: Rozk lad liczb otrzymanych metod

,

a Metropolisa dla rozk ladu Erlanga

20

Tablica 6: Ilo´s´

c odrzuconych liczb na 1000000 wylosowanych dla rozk ladu Erlanga

Metoda

losowanie 1

losowanie 2

losowanie 3

losowanie 4

´

srednia

von Neumanna

2992998

2978895

3000400

3002972

2993816.25

Metropolisa

51399

51296

50628

51275

51149.5

Kolejnym krokiem by lo sprawdzenie efektywno´sci metod, rozumianej jako stosunek liczby zaak-
ceptowanych zmiennych losowych do liczby wszystkich wygenerowanych zmiennych losowych.
Otrzymane wyniki dla 4 losowa´

n, a tak˙ze ´sredni stosunek dla obu metod zebrane zosta ly w

tabeli 7.

Tablica 7: Stosunek zaakceptowane/wylosowane dla rozk ladu Erlanga

Metoda

losowanie 1

losowanie 2

losowanie 3

losowanie 4

´

srednia

von Neumanna

25.0438392%

25.1326059%

24.9975002%

24.9814388%

25.038846025%

Metropolisa

95.1113707%

95.1206891%

95.1811678%

95.1225892%

95.1339542%

O ile w przypadku rozk ladu normalnego omawiannego w poprzednim podpunkcie r´

o˙znice w

ilo´sci generowanych liczb by ly du˙ze o tyle w tym wypadku s

,

a one wr

,

ecz olbrzymie. Przewaga

metody Metropolisa jest tu wyra´

znie widoczna. Co ciekawe, nie przek lada si

,

e ona bezpo´srednio

na szybko´s´

c dzia lania programu ani zu˙zycie pami

,

eci. Czas generacji przy u˙zyciu metody von

Neumanna jest o ok. 20-25% d lu˙zszy, a w wykorzystaniu pami

,

eci operacyjnej wogule nie wida´

c

r´

o˙znic mi

,

edzy tymi algorytmami.

Warto te˙z zauwa˙zy´

c, i˙z ilo´s´

c odrzuconych liczb w metodzie Metropolisa zale˙zna jest bezpo´srednio

od warto´sci parametru ∆. Aby otrzyma´

c podobn

,

a jak w poprzednim podpunkcie ilo´s´

c odrzu-

conych nale˙za lo by ustawi´

c ∆ ≈ 3.0, ale wtedy otrzymane wyniki r´

o˙zni

,

a si

,

e znacz

,

aco od

po˙z

,

adanego rozk ladu.

5.3

Dwuwymiarowy rozk lad normalny

Rozk ladem normalnym dwuwymiarowym nazywamy rozk lad o funkcji g

,

esto´sci prawdopodo-

bie´

nstwa postaci:

f(x

1

, x

2

) =

1

2π|Σ|

1/2

exp(−1/2(x

1

− µ

1

)Σ

−1

(x

1

− µ

1

))exp(−1/2(x

2

− µ

2

)Σ

−1

(x

2

− µ

2

))

gdzie: Σ jest macierz

,

a kowariancji zmiennych x

1

, x

2

, a |Σ| jest wyznacznikiem tej macierzy.

Macierz kowariancji definiujemy jako Σ =

Cov(x

1

, x

1

) Cov(x

1

, x

2

)

Cov(x

2

, x

1

)

Cov(x

2

x

2

)

!

Poniewa˙z Cov(x

1

, x

1

) = V ar(x

1

), a przyj

,

eto wariancj

,

e r´

own

,

a 1 (jak w rozk ladzie N(0, 1)) dla

obu zmiennych, macierz ta ma posta´

c Σ =

1 0
0 1

!

= Σ

−1

= I

2

czyli jest to macierz jednostkowa i |Σ| = 1.
Ostatecznie wi

,

ec g

,

esto´s´

c prawdopodobie´

nstwa dwuwymiarowego rozk ladu normalnego ( przyj-

muj

,

ac µ

1

= 0 i µ

2

= 0) ma posta´

c:

21

f (x

1

, x

2

) =

1

2π

exp(−

x

2
1

+ x

2
2

2

)

Inny, ciekawy spos´

ob otrzymania tego wyniku zaprezentowany jest w [17].

Tak jak i w poprzednich podpunktach wygenerowano 1000000 liczb, jednak tym razem ka˙zda
z tych liczb posiada dwie wsp´

o lrz

,

edne wi

,

ec by la to de facto generacja 2000000 liczb losowych.

Czas wykonywania metody Metropolisa wynosie l oko lo 13 sekund, a metody von Neumanna
oko lo 20. Jest to wzrost w stosunku do przypadku jednowymiarowego o oko lo 50%, czyli mniej
ni˙z si

,

e spodziewano. Oczywi´scie wp lyw na to mog ly mie´

c pewne zmiany w sposobie przydzie-

lania liczb do przedzia l´

ow, kt´

ore nieco przyspieszy ly prac

,

e programu.

Po otrzymaniu wynik´

ow wykonano wykresy przedstawiaj

,

ace rozk lad par liczb. Kolejno dla me-

tody von Neumanna (rys. 14) i Metropolisa (rys. 15). Wykre´slono tak˙ze teoretyczny rozk lad
g

,

esto´sci dla por´

ownania (rys. 16). Warto zauwa˙zy´

c, i˙z dla tego wykresu warto´s´

c funkcji zosta la

pomno˙zona przez i lo´s´

c liczb a nast

,

epnie przez obie skale podzia lu, czyli dwukrotnie przez 1:10.

Oczywi´scie wi

,

a˙ze si

,

e to z faktem, ˙ze nasz

,

a podstawow

,

a jednostk

,

a nie bedzie teraz odcinek jed-

nostkowy a kwadrat o boku 1. Same wykresy nie pokazuj

,

a zauwa˙zalnych odst

,

epstw rozk lad´

ow

otrzymanych od rozk ladu teoretycznego. Mo˙zna wi

,

ec uzna´

c i˙z obie metody s

,

a skuteczne w tym

wypadku.
Nie zamieszczono wykres´

ow otrzymanych liczb gdy˙z przy ich du˙zej liczbie wykres tr´

ojwymiarowy

stawa l sie zupe lnie nieczytelny i nie dawa l podstaw do ˙zadnych wniosk´

ow na temat otrzyma-

nych wynik´

ow. Oczywi´scie przedstawia l on naturalne rozszerzenie wykres´

ow 5 i 6 na przestrze´

n

tr´

ojwymiarow

,

a.

Dla metody von Neumanna stosowano parametr ∆ = 1.9 podobnie jak w podpunkcie pierw-
szym. Pozwala to na por´

ownanie otrzymanych wynik´

ow a jednocze´snie daje dobre efekty je´sli

chodzi o rozk lad otrzymanych liczb. Tabele 8 i 9 przedstawiaj

,

a odpowiednio ilo´s´

c odrzuconych

par liczb przy losowaniu 1000000 par oraz efektywno´s´

c obu metod rozumian

,

a jak w poprzednich

podpunktach (oczywi´scie tym razem dla par liczb).

Tablica 8: Ilo´s´

c odrzuconych par liczb na 1000000 wylosowanych par

Metoda

losowanie 1

losowanie 2

losowanie 3

losowanie 4

´

srednia

von Neumanna

9219146

9256037

9230236

9246219

9237909.5

Metropolisa

1074308

1069958

1071477

1070385

1071532

Jak wida´

c efektywno´s´

c obu metod jest zdecydowanie ni˙zsza ni˙z w obu poprzednich podpunk-

tach. Ma to oczywi´scie zwi

,

azek z faktem losowania pary liczb zamiast jednej liczby. Mo˙zna

to wyt lumaczy´

c jako zmniejszenie si

,

e stosunku dost

,

epnej dla liczb cz

,

e´sci przestrzeni wykresu

z oko lo

1
3

w podpunkcie pierwszym do oko lo

1
9

w tym przypadku. Praktycznie odpowiada

to proporcjom pomi

,

edzy opdowiednimi warto´sciami dla rozk ladu normalnego jednowymiaro-

wego i dwuwymiarowego dla metody von Neumanna. Wystarczy por´

owna´

c pierwsze wiersze z

tabel 3 i 9 aby zauwa˙zy´

c, ˙ze w przypadku dwuwymirowym mamy oko lo 3 razy mniejsz

,

a efek-

tywno´s´

c. Z tego samego por´

ownania wida´

c te˙z, ˙ze spostrze˙zenie to nie odnosi si

,

e do metody

Metropolisa. Wynika to z faktu, i˙z metoda ta opiera si

,

e na zasadzie b l

,

adzeniea przypadko-

wego i istnieje bardzo ma le prawdopodobie´

nstwo zaakceptowania krokupoza g l´

ownym obsza-

rem ´s lupka”widocznego na wykresie 15. Tym samym praktycznie nigdy nie otrzymujemy liczb
spoza zakresu (−3, 3) na obu osiach gdy˙z wykonie dw´

och kolejnych krok´

ow w stron

,

e brzegu

obszaru wykresu jest bardzo nieprawdopodobne. Dlatego te˙z wi

,

ekszo´s´

c krok´

ow wykonywana

jest w stron

,

e s

,

asiedzwa ponktu (0, 0).

Wa˙znym wnioskiem jaki wyci

,

agn

,

a´

c mo˙zna z przedstawionych tabel jest, ˙ze przy wi

,

ekszej liczbie

wymiar´

ow metoda von Neumanna staje si

,

e bardzo czaso- i zasoboch lonna. Natomiast metoda

22

Rysunek 14: Rozk lad liczb otrzymanych metod

,

a von Neumanna dla rozk ladu normalnego dwu-

wymiarowego

Rysunek 15: Rozk lad liczb otrzymanych metod

,

a Metropolisa dla rozk ladu rozk ladu normalnego

dwuwymiarowego

23

Rysunek 16: Rozk lad normalny dwuwymiarowy

Metropolisa z odpowiednio dobranym parametrem ∆ daj wyniki r´

ownie dobre przy znacznie

mniejszym zu˙zyciu zasob´

ow maszynowych.

Tablica 9: Stosunek zaakceptowane/wylosowane dla rozk ladu normalnego dwuwymiarowego

Metoda

losowanie 1

losowanie 2

losowanie 3

losowanie 4

´

srednia

von Neumanna

9.7855535%

9.7503548%

9.7749456%

9.7596977%

9.7676379%

Metropolisa

48.2088484%

48.310159%

48.2747334%

48.3001954%

48.27984955%

24

6

Podsumowanie

Por´

ownania obu metod dokonano dla trzech r´

o˙znych rozk lad´

ow:

1. standardowego rozk ladu normalnego N(0,1)

2. rozk ladu Erlanga

3. rozk ladu normalnego dwuwymiarowego (binormal )

Otrzymano dla ka˙zdego z nich zadowalaj

,

ace wyniki je´sli chodzi o zgodno´s´

c z rozk ladem teo-

retycznym. Testy przeprowadzono tylko dla pierwszego przypadku, da ly one jednak wyniki
pozytywne wi

,

ec mo˙zna przypuszcza´

c, i˙z w pozosta lych przypadkach by loby podobnie. Wa˙znym

wynikiem otrzymanym dla ka˙zdego z przypadk´

ow by la ich efektywno´s´

c czyli stosunek liczby

zaakceptowanych zmiennych losowych do liczby wszystkich wygenerowanych zmiennych loso-
wych (dla przypadku trzeciego jako zmienn

,

a traktujemy punkt czyli par

,

e liczb). Por´

ownanie

efektywno´sci metod w zale˙zno´sci od generowanego rozk ladu zawiera tabela 10.

Tablica 10: Por´

ownanie efektywno´sci metod

Metoda

rozk lad N(0, 1)

rozk lad Erlanga

rozk lad normalny 2D

´

srednia

von Neumanna

31.262156525%

25.038846025%

9.7676379%

22.02288015%

Metropolisa

73.92361635%

95.1339542%

48.27984955%

72.4458067%

Trzeba pami

,

eta´

c, ˙ze dla rozk ladu Erlanga warto´s´

c parametru ∆ = 0.5 dla metody Metropolisa,

czyli mniejsza ni˙z w pozosta lych przypadkach co spowodowa lo mniejsz

,

a ilo´s´

c odrzuconych liczb.

Dla warto´sci ∆ = 1.9, czyli takiej jak w przypakdu 1. i 3. efeektywno´s´

c wynios laby oko lo 81%

a ´srednia dla taj metody spad laby do ok 67%. Co ciekawe badano r´

ownie˙z parametr ∆ = 4.0

dla 3. przypadku co powinno odpowiada´

c takiej samej dowolno´sci losowania jak w metodzie

von Neumanna ale ilo´s´

c odrzuconych liczb by la i tak o oko lo 40% mniejsza ni˙z w tej metodzie.

Og´

olnie jednak trzeba oceni´

c oba algorytmy jako w miar

,

e szybkie oraz, co wa˙zne, latwe w

implementacji co czyni je bardzo dobrymi generatorami liczb losowych o zadanym rozk ladzie.
Metoda Metropolisa jest efektywniejsz

,

a jednak wymaga dobrego wyboru parametru ∆ co dla

du˙zych przedzia l´

ow nie jest zadaniem latwym.

25

Spis rysunk´

ow

1

Metoda von Neumanna . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

2

Generator mieszany . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

13

3

Funkcja rand()

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

13

4

Marsenne Twister . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

13

5

Liczby otrzymane metod

,

a von Neumanna dla rozk ladu N(0,1) w 2D . . . . . . .

15

6

Liczby otrzymane metod

,

a Metropolisa dla rozk ladu N(0,1) w 2D . . . . . . . . .

15

7

Rozk lad liczb otrzymanych metod

,

a von Neumanna dla g(x) = 1/2 * sqrt(2e/π)

* exp(-|x|) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

16

8

Rozk lad liczb otrzymanych metod

,

a von Neumanna dla rozk ladu N(0,1) w 2D . .

17

9

Rozk lad liczb otrzymanych metod

,

a Metropolisa dla rozk ladu N(0,1) w 2D . . . .

17

10

Liczby otrzymane metod

,

a von Neumanna dla rozk ladu Erlanga

. . . . . . . . .

19

11

Liczby otrzymane metod

,

a Metropolisa dla rozk ladu Erlanga . . . . . . . . . . .

19

12

Rozk lad liczb otrzymanych metod

,

a von Neumanna dla rozk ladu Erlanga

. . . .

20

13

Rozk lad liczb otrzymanych metod

,

a Metropolisa dla rozk ladu Erlanga . . . . . .

20

14

Rozk lad liczb otrzymanych metod

,

a von Neumanna dla rozk ladu normalnego

dwuwymiarowego . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

23

15

Rozk lad liczb otrzymanych metod

,

a Metropolisa dla rozk ladu rozk ladu normal-

nego dwuwymiarowego . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

23

16

Rozk lad normalny dwuwymiarowy . . . . . . . . . . . . . . . . . . . . . . . . . .

24

26

Spis tablic

1

´

Srednia liczb na przedziale [0,1) . . . . . . . . . . . . . . . . . . . . . . . . . . .

12

2

Ilo´s´

c odrzuconych liczb na 1000000 wylosowanych dla N(0,1) w 2D . . . . . . . .

14

3

Stosunek zaakceptowane/wylosowane dla N(0,1) w 2D . . . . . . . . . . . . . . .

14

4

Wyniki testu chi-kwadrat dla obu metod . . . . . . . . . . . . . . . . . . . . . .

16

5

Wyniki testu Ko lmogorowa-Smirnowa dla obu metod . . . . . . . . . . . . . . .

18

6

Ilo´s´

c odrzuconych liczb na 1000000 wylosowanych dla rozk ladu Erlanga . . . . .

21

7

Stosunek zaakceptowane/wylosowane dla rozk ladu Erlanga . . . . . . . . . . . .

21

8

Ilo´s´

c odrzuconych par liczb na 1000000 wylosowanych par . . . . . . . . . . . . .

22

9

Stosunek zaakceptowane/wylosowane dla rozk ladu normalnego dwuwymiarowego

24

10

Por´

ownanie efektywno´sci metod . . . . . . . . . . . . . . . . . . . . . . . . . . .

25

27

Literatura

[1] Ryszard Zieli´

nski: ”Generatory liczb losowych: Programowanie i testowanie na maszynach

cyfrowych” WNT, Warszawa 1979

[2] http://www.wikipedia.org/

[3] J. von Neumann: ”Various techniques used in connection with random digits. Monte Carlo

methods” Nat. Bureau Standards, 12 (1951), pp. 36 - 38

[4] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller:

,

Equation

of state calculation by fast computing machines.” Journal of Chemical Physics, 21(6):1087
- 1092, 1953

[5] Wyk lad Metody Obliczeniowe fizyki prof. dr hab. Janusz Adamowski

[6] http://www.zftik.agh.edu.pl/mof/lab/MOF B 4 2.pdf

[7] http://fatcat.ftj.agh.edu.pl/yaq/random.pdf

[8] http://sphere.pl/nadolski/kms/w/w3.pdf

[9] http://omega.albany.edu:8008/cdocs/summer99/lecture4/l4.ps

[10] http://www.cplusplus.com/

[11] http://www.gnuplot.info/

[12] http://www.azillionmonkeys.com/qed/random.html

[13] http://www.math.sci.hiroshima-u.ac.jp/ m-mat/MT/emt.html

[14] http://www.zftik.agh.edu.pl/mof/lab/a3.pdf

[15] http://www.agner.org/random/randomc.htm

[16] http://home.agh.edu.pl/ bartus/

[17] http://mathworld.wolfram.com/BivariateNormalDistribution.html

28