Generatory liczb pseudolosowych

Generatory liczb pseudolosowych

Plan wykładu:

1.Generatory o rozkładzie równomiernym a) liniowe

b) kombinowane – generator uniwersalny c) nieliniowe

2. Generatory o dowolnym rozkładzie prawdopodobieństwa a) metoda odwracania dystrybuanty

b) metoda eliminacji

c) superpozycja rozkładów d) rozkład dyskretny

3.Generatory o rozkładach wielowymiarowych a) Rozkład równomierny na sferze i kuli w R^M b) Dwuwymiarowy rozkład normalny

4.Testowanie generatorów

a) Testy zgodności z rozkładem: ^, OPSO

b) Testy zgodności z rozkładem statystyk: pozycyjne, test sum c) Testy kombinatoryczne

2 Generatory o rozkładzie równomiernym

Najprostsze generatory liczb losowych to generatory fizyczne wykorzystujące:

1) szumy układów elektronicznych 2) Promieniotwórczość

Zalety: dostajemy ciągi liczb losowych (niezależne, nieskorelowane)

Wady: wymagana ciągła kalibracja (testowanie parametrów), kłopoty techniczne z obsługą, brak powtarzalności serii.

Własności generatorów komputerowych a) łatwość obsługi

b) możliwość generowania dowolnego rozkładu c) dowolna liczba wymiarów

d) Powtarzalność ciągów generowanych liczb generatory

liczb losowych

fizyczne programowe

(komputerowe)

Jak pracują komputerowe generatory liczby?

1) Tworzony jest ciąg liczb nieujemnych (naturalnych)

o rozkładzie równomiernym według wybranego algorytmu. Liczby ograniczone są od góry przez reprezentację, np. dla k=32 bitowej reprezentacji liczby całkowitej bez znaku, kres górny

T_a – okres aperiodyczności ciągu T₀ – okres ciągu

2) Aby uzyskać liczby „rzeczywiste” dokonujemy przekształcenia

3) Dokonujemy kolejnej transformacji ciągu aby uzyskać ciąg o zadanym rozkładzie

prawdopodobieństwa (normalny, wielomianowy, etc.)

X

; X

; X

; : : : ; X

m = 2

U = X

m ) U

2 (0; 1]

X

; : : : ; X

| {z }

Ta

X

; X

; : : : ; X

| {z }

To

; X

; : : :

3 Generatory liniowe

Generatory liniowe tworzą ciąg liczb według schematu:

gdzie:

a₁, a₂,...,a_k, c, m – parametry generatora (ustalone liczby)

Operację

nazywamy dzieleniem modulo a jej wynikiem jest reszta z dzielenia liczb całkowitych a i b.

Lub inaczej: r jest kongruentne do A modulo n jeśli n jest dzielnikiem a-r.

Generatory wykorzystujące operację dzielenia modulo to generatory kongruentne lub

kongruencyjne.

Przykład

19 mod 6 =1 15 mod 6 =3 18 mod 6 =0 14 mod 6 =2 17 mod 6 =5 13 mod 6 =1 16 mod 6 =4 12 mod 6 =0

X

= (a

X

+ a

X

+ : : : +a

X

+ c) mod m

r = (a mod n); a; n; r 2 Z

a ´ r mod n ) r = a ¡ j a n

k n

Aby wygenerować ciąg liczb pseudolosowych należy zdefiniować jego parametry.

Liczby

nazywamy ziarnem generatora (seed). Dla bardziej rozbudowanych generatorów liczby te otrzymujemy z innego generatora lub np. używając zegara

systemowego (X₀).

Najprostszy generator liniowy ma dwie odmiany a) generator multiplikatywny gdy

b) generator mieszany gdy

Maksymalny okres generatora liniowego to (m-1) Generator multiplikatywny

c 6= 0 c = 0

X

; X

; X

; : : : ; X

X

= aX

mod m k

=

¹ aX

m

º

; i ¸ 1 X

= aX

¡ mk

X

= a

X

¡ mk

¡ mk

a

X

= a

X

¡ mk

¡ mk

a ¡ mk

a

: : : : : : : : : : : :

X

= a

X

¡ m(k

+ k

a + : : : + k

a

)

4 Ostatnie równanie można zapisać w postaci

skąd wynika, że wybór X₀ determinuje wszystkie liczby w generowanym ciągu (a i m są ustalone) – uzyskany ciąg liczb jest deterministyczny

Przykład.

Generator multiplikatywny

X

= a

X

mod m

Xi

a=1 2 3 4 5 6 7 8 9 10

i=0 1 1 1 1 1 1 1 1 1 1

1 1 2 3 4 5 6 7 8 9 10

2 4 9 5 3 3 5 9 4 1

3 8 5 9 4 7 2 6 3

4 5 4 3 9 9 3 4 5

5 10 1 1 1 10 10 10 1

6 9 5 4 3

7 7 8 6 2

8 3 4 9 5

9 6 2 8 7

10 1 1 1 1

X

= aX

mod 11; X

= 1

5 Okres generatora multiplikatywnego

Maksymalny okres generatora multiplikatywnego wynosi

Gdy m jest liczbą pierwszą a p jest czynnikiem pierwszym liczby (m-1).

Przykład

Wykorzystujemy liczby Marsenne'a (które dość często są liczbami pierwszymi)

Okres generatora

Liczby

występują dokładnie 1 raz w pojedynczym okresie generatora.

Odległość pomiędzy najbliższymi sąsiadami

T = min fi : X

= X

; i > 0 g

a

6= 1 mod m

m = 2

¡ 1

a ´ 1 mod 4

a ´ 1 mod p

a = 69069 c = 1 m = 2

W przypadku generatora mieszanego maksymalny okres uzyskamy po spełnieniu warunków:

a) liczby c i m nie mają wspólnych dzielników b) jeżeli 4 jest dzielnikiem m to

c) dla każdego dzielnika pierwszego p liczby m

Przykład.

T = 2

¡ 2

p = 31 ) m = 2

; a = 7

U = fU

; 1 · i · T g

1

2

¡ 1 = 4:657 £ 10

6 Rozkład przestrzenny ciągu

Wadą generatorów multiplikatywnych jest

nierównomierne pokrycie d-wymiarowej kostki (I^d).

Generowane liczby lokalizują się na

hiperpłaszczyznach, których położenie uzależnione jest od parametrów generatora.

Przykład.

X

= aX

mod 11

a) X₀=1, a=2 b) X₀=1, a=8

(U

; U

; : : : ; U

); (U

; U

; : : : ; U

); : : :

(U

; U

; : : : ; U

); (U

; U

; : : : ; U

); : : :

7 Parametry statystyczne generatora o rozkładzie

równomiernym w (0,1)

Jeśli generowany ciag liczb jest niezależny to wartość oczekiwana (średnia) powinna wynosić

natomiast wariancja jest równa

Ponadto współczynniki autokorelacji elementów ciągu powinny wynosić 0.

Jeśli parametry statystyczne generatora (ciągu generowanych przez niego liczb) odbiegają od powyższych wartości to jest on nieprzydatny (lub warunkowo przydatny).

¹ = Z

0

xdx = x

2 j

= 1 2

¹

¹ = 1 N

X

x

¾

= Z

0

(x ¡ ¹)

dx = 1 12

¹

¾ = 1 N

X

(x

¡ ¹¹)

Przykład.

Przeanalizujmy parametry statystyczne L bitowego generatora

Można nim wygenerować 4 ciągi liczb o okresie

Jeden z nich:

Ciąg liczb pseudolosowych

jest permutacją liczb

Liczby znormalizowane

Średnia

Wariancja

X

= aX

mod 2

c ´ 3 mod 8; X

´ 1; 3; 9; 11 mod 16 X

; X

; X

; : : :

8j + 1; 8j + 3; j = 0; 1; : : : ; 2

¡ 1 T = 2

U

= X

2

¾

= 1

12 ¡ 13 3 ¢ 2

¹

¹ = 1

2 ¡ 1

2

8 Funkcja autokorelacji opisuje zależność

elementów ciągu od wyrazów poprzednich.

Definicja

Oraz wzór dla ciągu skończonego

Inaczej: opisuje związek pomiędzy elementami dówch szeregów – danego i przesuniętego o r.

Efektywnie funkcję autokorelacji można badać przy użyciu FFT – splot dwóch wektorów.

Dla analizowanego generatora mieszanego oraz dla r=1 współczynnk autokorealcji można oszacować z poniższej relacji

A = 1

a ¡ 6c

am (1 ¡ c m )

B = a

m

R

2 [A ¡ B; A + B]

Przykłady generatorów liniowych

Ich okresy są maksymalne tj. równe (m-1) Generatory na rejestrach przesuwnych Definiujemy ciąg bitów b_i otrzymywanych rekurencyjnie

gdzie:

są stałymi binarnymi , a stałe b_i

tworzą ciąg inicjujący.

X

= (1176X

+ 1476X

+ 1776X

) mod (2

¡ 5) X

= 2

(X

+ X

+ X

) mod (2

¡ 5)

X

= (1995X

+ 1998X

+ 2001X

) mod (2

¡ 849) X

= 2

(X

+ X

+ X

) mod (2

¡ 1629)

b

= (a

b

+ : : : + a

b

) mod 2 i = k + 1; k + 2; : : :

a

; a

; : : : ; a

2 f0; 1g

b

; b

; : : : ; b

2 f0; 1g R ~

= lim

n!1

1 N

N

X

X

X

R

= 1

(N ¡ r)¾

N

X

(X

¡ ¹)(X

¡ ¹)

9 Inny sposób zapisu relacji rekurencyjnej

wykorzystuje operator xor

Który można zdefiniować

Jeśli założymy

to relację rekurencyjną można zapisać przy użyciu xor

Jakie są własności takiego ciągu bitów?

Ciąg jest okresowy o okresie

nieprzekraczającym 2^k a dokładniej (2^k-1) ze względu na wyrzucenie układu bitów złożonych z samych zer.

Praktyczną realizacją tego pomysłu jest

algorytm wykorzystujący tylko dwa elementy ciągu:

a b a xor b

0 0 0

0 1 1

1 0 1

1 1 0

a xor b = (a + b) mod 2 a

= a

= : : : = a

= 1

b

= b

xor b

xor : : : xor b

b

= b

xor b

; p > q; p; q 2 N

Wykorzystujemy ciąg bitów do obliczenia liczby pseudolosowej z przedziału (0,1]

(generator Tauswortha)

gdzie: s jest ustaloną całkowitą liczbą nieujemną a) jeśli s<L to do utworzenia U_i oraz U_i+1

wykorzystywane są elementy tego samego podciągu

b) jeśli s=L to U_i oraz U_i+1 są tworzone z rozłącznych fragmentów ciągu globalnego

Ciąg bitów łatwo generuje się przy użyciu rejestrów przesuwnych oraz bramek logicznych (xor) – łatwa implementacja w języku C.

U

=

X

2

b

= 0:b

: : : b

i = 0; 1; 2; : : :

10 Generator Fibonacciego

Punktem wyjścia do konstrukcji generatora są liczby Fibonacciego

a dokładniej ciąg reszt tego ciągu

Powyższy ciąg reszt ma rozkład równomierny ale nie spełnia testów niezależności.

Jego modyfikacja

nie posiada już tej wady.

Okres generatora zależy od operacji. Dla

f

= f

= 1

f

= f

+ f

¦ = +; ¡; ¤; xor

T = 8 >

> <

> >

:

(2

¡ 1)2

F (r; s; +) (2

¡ 1)2

F (r; s; ¡) (2

¡ 1)2

F (r; s; ¤)

(2

¡ 1) F (r; s; xor) m = 2

Generator charakteryzuje się dużym okresem ale jest wolniejszy np. względem generatorów

multiplikatywnych co obecnie nie jest już dużą wadą.

Kombinacje generatorów

Zakładamy że dysponujemy zmiennymi losowymi X i Y określonymi na zbiorze S={1,2,...,n}

z rozkładami prawdopodobieństwa

Z tych liczb możemy utworzyć wektory

Za normę wektora przyjmiemy p-normę

Miarą „odległości” danego rozkładu od rozkładu rónomiernego będzie wówczas wyrażenie

qqq = (q

; q

; : : : ; q

)

krk =

Ã

X

i=1

r

!

rrr = (r

; r

; : : : ; r

) P fX = ig = r

P fY = ig = q

i = 1; 2; : : : ; n

±(X) = k(r

; r

; : : : ; r

) ¡ (1=n; 1=n; : : : ; 1=n)k X

= X

+ X

mod m; i ¸ 2

X

= X

¦ X

mod m

i ¸ r; r > s ¸ 1

11 Dla określonego działania na zbiorze S tj.

rozkład nowej zmiennej losowej

będzie bliższy rozkładowi równomiernemu niż rozkłady zmiennych X i Y

Nowy ciąg ma lepsze własności statystyczne a także większy okres. Jeżeli ciąg

ma okres T₁, a ciąg

ma okres T₂ i są liczbami pierwszymi to wtedy nowy ciąg liczb

ma okres równy T₁T₂.

Generator uniwersalny

Daje jednakowe wyniki na dowolnym komputerze. Jego działanie oparte jest na wykorzystaniu kombinacji kilku generatorów.

+; ¡; ¤; xor X ¦ Y

±(X ¦ Y ) · minf±(X); ±(Y )g

X

; X

; : : : ; X

Y

; Y

; : : : ; Y

X

¦ Y

Pierwszy z nich jest generatorem Fibonacciego

działanie jest zdefiniowane następująco

Inicjalizację generatora tj.wyznaczenie ciągu

przeprowadzamy przy pomocy ciągu bitów (24-bitowa mantysa, 16 bitowe liczby całkowite) tj.

Ciąg bitów generujemy przy użyciu dwóch generatorów

V

= V

¦ V

V

2 [0; 1)

x ¦ y = x ¡ y; x ¸ y F (97; 33; ¦)

x ¦ y = x ¡ y + 1; x < y V

; V

; : : : ; V

V

= 0:b

b

: : : b

V

= 0:b

b

: : : b

y

= (y

y

y

) mod 179 z

= (52z

+ 1) mod 169 b

=

½ 0; gdy y

z

mod 64 < 32

b

= 1; w pozostalych przypadkach

U

= V

¦ c

12 Generatory muszą być zainicjowane

Okres generatora kombinowanego: 2¹²⁰.

Drugi generator (c_n) o rozkładzie równomiernym w (0,1) jest zdefiniowany następująco

Okres tego generatora to 2²⁴-3.

Okres generatora uniwersalnego (U_n): 2¹⁴⁴

y

; y

; y

2 f1; 2; : : : ; 178g z

2 f0; 1; : : : ; 168g

c

= c

¦

µ 7654321 16777216

¶

n ¸ 2

c

= 362436 16777216 c ¦ d = c ¡ d; c ¸ d c ¦ d = c ¡ d + 362436

16777216 ; c < d c; d 2 [0; 1)

Generator ULTRA – oparty na odejmowaniu z pożyczką

Wykorzystujemy w nim operację

której wynikiem jest

Parametr c jest tzw. bitem przeniesienia. Jego początkowa wartośc wynosi 0.

W generatorze ULTRA przyjęto parametryzację

Elementy ciągu wyznaczamy rekurencyjnie

Generator inicjalizujemy L-bitowymi liczbami

x ª y mod m

½ x ¡ y ¡ c + m x ¡ y ¡ c < 0 ! c = 1 x ¡ y ¡ c x ¡ y ¡ c ¸ 0 ! c = 0

c 2 f0; 1g

m = 2

L ¸ 32; r; s; r > s

X

= X

ª X

mod m X

= b

b

: : : b

X

= b

b

: : : b

X

= : : :

13 Ciąg bitów

tworzymy ze znaków liczb

które generujemy następująco:

1) Podajemy dwie liczby:

2) Kolejno obliczamy

gdzie: k₁, k₂ są ustalonymi liczbami.

Operacja

oznacza przesunięcie bitów w v_i w prawo i wypełnienie najstarszych k₁ bitów zerami.

b

; b

; : : :

w

; w

; : : :

u

; v

2 (0; m)

u

= ¸u

mod m v

= (v

>> k

) xor v

v

= (v

<< k

) xor v

w

= u

xor v

v

>> k

Jeśli

jest liczbą pierwszą to okres generatora jest równy najmniejszej liczbie naturalnej k dla której spełnione jest poniższe równanie

W generatorze ULTRA dla parametrów

Okres generatora wynosi 10³⁴⁶.

Generator jest bardzo szybki. Pozwala generować liczby całkowite: 32, 31, 8, 7 i 1 bitowe.

M = m

¡ m

+ 1

m

mod M = 1

L = 32 r = 37 s = 24

¸ = 69069

14 Generatory nieliniowe

Zaletą generatorów nieliniowych jest to że ciągi generowanych przez nie liczb nie układają się na hiperpłaszczyznach tak jak w przypadku

generatorów liniowych.

Generator nielinowy Eichenauera-Lehna

gdzie: m jest liczbą pierwszą

Odwrotność modulo liczymy następująco

Uzyskujemy w ten sposób ciąg liczb

o rozkładzie równomiernym.

c = 0 ) c

mod m = 0

c 6= 0 ) c ¢ (c

mod m) = 1 c

= c

mod m

X

2 f0; 1; : : : ; m ¡ 1g

U

= X

m 2 [0; 1)

Generator Eichenauera-Hermanna

W generatorze tym kolejne elementy ciągu są niezależne od poprzednich – cecha szczególnie przydatna w obliczeniach równoległych.

Dla

generator osiąga okres T=m.

Aby generator osiągał maksymalny okres równy m musi być spełniony jeszcze warunek

jest najmniejszą liczbą całkowitą dla której zachodzi

a 2 f1; 2; : : : ; mg

m

¡ 1

z

= 1 mod (z

¡ bz ¡ a) X

= (aX

+ b) mod m; i = 0; 1; : : :

X

= (a(i + i

) + b)

mod m; i = 0; 1; : : :

15 Generatory o dowolnym rozkładzie

prawdopodobieństwa

Metoda odwracania dystrybuanty Dystrybuanta określonego rozkładu prawdopodobieństwa jest funkcją

niemalejącą i prawostronnie ciągłą

Dystrybuanta jednoznacznie definiuje rozkład prawdopodobieństwa.

Związek pomiędzy dystrybuantą a gęstością prawdopodobieństwa f(x):

Przykład.

Rozkład wykładniczy

F : R ! R

x!¡1

lim F (x) = 0

x

lim

F (x) = 1

F (x) = Z

¡1

f (t)dt

F (x) =

½ 1 ¡ e

x ¸ 0 0 x < 0

Rozkład normalny

F (x) = Z

¡1

1

¾ p

2¼ exp µ

¡ (t ¡ ¹)

2¾

¶ dt

Rys. Dystrybuanta rozkładu normalnego.

f (x) = e

16 U jest zmienną losową o rozkładzie

równomiernym w przedziale (0,1).

Nową zmienną losową będzie

Generujemy więc ciąg liczb pseudolosowych

który przekształcamy w ciąg

Liczby X_i mają rozkład prawdopodobeństwa o dystrybuancie F.

Odwracanie dystrybuanty można wykorzystać także w przypadku rozkładów dyskretnych.

Np. ciąg zmiennych

o rozkładzie

X = F

(U )

P fX · xg = P fF

(U ) · xg

= P fU · F (x)g

= F (x)

U

; U

; : : : ; U

2 (0; 1)

X

; X

; : : : ; X

2 (¡1; 1)

X

; X

; : : : ; X

p

= P fX = kg; k = 0; 1; 2; : : :

Można wygenerować przy użyciu ciągu

korzystając ze wzoru

Odwracanie dystrybuanty sprawia duże trudności numeryczne.

Przykład.

Odwracanie dystrybuanty dla rozkładu normalnego

gdzie: L(t) i M(t) są wielomianami stopnia 4 zmiennej t.

Zaletą tej metody jest to, że do wygenerowania jednej nowej zmiennej losowej X_n potrzeba tylko jednej zmiennej pierwotnej U_n.

U

; U

; : : : ; U

2 (0; 1)

X

= min (

k : U

· X

i=0

p

)

; n = 1; 2; : : :

©

(u) =

½ g(u) 10

< u < 0:5

¡g(1 ¡ u) 0:5 · u · 1 ¡ 10

g(u) = t ¡ L(t)

M (t) ; t = p

¡2ln u

17 Metoda eliminacji (von Neumann)

Chcemy wygenerować ciąg zmiennych losowych o gęstości prawdopodobieństwa f w przedziale [a,b].

Wartość f jest w przedziale [a,b] ograniczona od góry przez stałą d.

Sposób otrzymania ciągu zmiennych losowych o rozkładzie f(x) jest następujący:

1) Losujemy dwie zmienne o rozkładzie równomiernym

2) jeżeli

3) gdy powyższy warunek nie jest spełniony wówczas odrzucamy parę U₁, U₂

4) wykonujemy czynności 1-3 aż do uzyskania odpowiednio licznego ciągu

Wygenerowana zmienna losowa X ma rozkład prawdopodobieństwa f.

Oznaczenia

- pole powierzchnii zbioru C (miara Lebesgue'a): l₂(C) - zbór A

Ponieważ f(x) jest gęstością prawdopodobieństwa, więc

U

2 [a; b] U

2 [0; d]

U

· f(U

) ) X = U

A = f(x; u) : a · x · b; 0 · u · f(x)g

l

(A) = 1

18 Pole powierzchnii prostokąta (a,b)x(0,d)

Przy tych założeniach możemy stwierdzić że dowolnie wybrany punkt  ma rozkład

równomierny na zbiorze

jeżeli dla podzbiorów C tego zbioru zachodzi

Tw. Jeśli (X₁,U₁), (X₂,U₂),... jest ciągiem punktów losowych o rozkładzie

równomiernym w prostokącie (a,b)x(0,d) oraz (X_i,U_i) będzie pierwszym punktem

wpadajacym do zbioru A to punkt ten ma rozkład równomierny w A.

Dowód. Rozważamy podzbiór B zbioru A

l

((a; b) £ (0; d)) = (b ¡ a)d

D ½ R

P f! 2 Cg = l

(C) l

(D)

P f(X; U) 2 Bg = X

i=1

P f(X

; U

) = 2 A; : : : ; (X

; U

) = 2 A; (X

; U

) 2 Bg

= X

i=1

µ

1 ¡ l

(A) (b ¡ a)d

¶

l

(B)

(b ¡ a)d = l

(B) l

(A)

Tw. Jeżeli U ma rozkład równomierny w [0,1] a X ma rozkład o gęstości f(x) oraz X i U są niezależne, to punkt losowy (X,Uf(X)) ma rozkład równomierny w zbiorze A.

Powyższe rozważania można uogólnić na większą liczbę wymiarów.

Ponadto zamiast stałej d ograniczającej f(x) od góry można wprowadzić funkcję ograniczającą cg(x), wówczas warunkiem akceptacji jest

gdzie: c jest stałą, a g(x) rozkładem zbliżonym do f(x).

Zabieg taki opłaca się wykonać wówczas gdy

wyznaczenie ciągu liczb o gęstości g(x) jest łatwe i szybkie. Funkcję g(x) nazywamy wtedy rozkładem dominującym.

Jak dobrać stałą c?

cg(X)U · f(X); X 2 [a; b]; U 2 [0; 1]

c = sup

x

f (x) g(x) P fU cg(X) · f (X)g = Z

R^k

g(x)d

x

Z

du = f (x)

cg(x)

19 Superpozycja rozkładów

Naszym zadaniem jest uzyskanie ciągu liczb pseudolosowych o gęstości f(x), którą możemy wyrazić w postaci

gdzie: g_t(x) oraz h(t ) są również gęstościami prawdopodobieństwa – funkcje znane.

Rozkład prawdopodobieństwa f(x) nazywamy rozkładem złożonym.

Algorytm generowania liczb o rozkładzie f(x) jest następujący:

1) Generujemy zmienną T o rozkładzie h

2)Dla wartości t zmiennej losowej T generujemy zmienną losową X dla rozkładu gęstości g_t(x) W praktyce całkę zastępuje się sumą

f (x) =

Z

¡1

g

(x)h(t)dt

f (x) = X

i=1

p

g

(x) = X

i=1

p

g

(x)

Przedział [a,b] w którym generujemy ciąg zmiennych z rozkładem f(x) dzielimy na sumę K rozłącznych podprzedziałów. W każdym z nich (A_i) wyznaczamy

oraz

gdzie: 1_Ai jest funkcją przynależności do podzbioru A_i. Algorytm generacji ciagu zmiennej X jest wówczas następujący:

1) Losujemy zmienną losową

2) Dla wygenerowanej wartości i zmiennej I generujemy X z rozkładem gęstości

Dla dostatecznie wąskich przedziałów można h(t) przybliżyć wielomianem niskiego stopnia.

p

= Z

A_i

f (x)dx

g

(x) = 111

f (x) p

I 2 f1; 2; : : : ; Kg

g

(x) = 111

f (x) p

p

¸ 0;

X

p

= 1

20 Rozkład dyskretny

a) metoda odwracania dystrybuanty W rozkładzie dyskretnym mamy określone prawdopodobieństwo wylosowania danej liczby

wraz z warunkiem

Algorytm generowania ciągu zmiennych o rozkładzie dyskretnym:

1) X=0, S=p₀

2) Losujemy zmienną U o rozkładzie równomiernym z przedziału [0,1]

3) Sprawdzamy warunek

a)Iteracyjnie obliczamy

dopóki warunek jest spełniony

b) W przeciwym wypadku akceptujemy X 4) Kroki 1-3 wykonujemy aż do uzyskania odpowiednio licznego ciągu X₁,X₂,X₃,...

P fX = kg = p

; k = 0; 1; 2; : : : ; M

X

p

= 1

U > S

X = X + 1; S = S + p

b) Metoda równomiernego rozbicia przedziału Zakładamy że zmienna losowa X przyjmuje określone wartości z pewnym prawdopodobieństwem

Przedział (0,1) dzielimy na K+1 podprzedziałów o jednakowej długości

Zmienna U wpada do przedziału

Konstruujemy dwa ciągi

µ i ¡ 1

K + 1 ; i K + 1

¶

; i = 1; 2; : : : ; K + 1

[(K + 1)U + 1]

P fX = kg = p

; k = 0; 1; 2; : : : ; K

q

= 0

g

= max

½

j : q

< i K + 1

¾

; i = 1; 2; : : : ; K + 1 q

=

X

p

; i = 0; 1; : : : ; K

21 Algorytm

1) Generujemy zmienną U o rozkładzie równomiernym w (0,1)

2) Obliczamy

3) Iteracyjnie obliczamy

dopóki jest spełniony warunek

4) Jeśli

to akceptujemy X

5) Kroki 1-4 powtarzamy aż do uzyskania odpowiednio licznego ciągu X₁,X₂,...

Powyższy algorytm zapewnia to, że warunek

będzie sprawdzany conajwyżej dwukrotnie.

X = [(K + 1)U + 1]

X = X ¡ 1 q

> U X = g

+ 1

q

< U

q

> U

Generatory o rozkładach wielowymiarowych Zadanie można sformułować następująco:

Należy wygenerować ciąg wielowymiarowych zmiennych losowych

których rozkład prawdopodobieństwa ma gęstość

Do generacji takiego ciągu można stosować metodę eliminacji czy superpozycji rozkładów. Przy użyciu metody elminacji w najprostszej postaci pojawiają się problemy.

Przykład.

Określić prawdopodobieństwo akceptacji

wielowymiarowej zmiennej losowej o rozkładzie równomiernym na kuli jednostkowej (K_m(0,1)) . Algorytm.

Losujemy m zmiennych niezależnych o rozkładzie równomiernym w (-1,1) i konstruujemy zmienną wielowymiarową

Zmienną akceptujemy jeśli

X

X X = (X

; X

; : : : ; X

)

f (x

; x

; : : : ; x

)

U U U = (U

; U

; : : : ; U

)

kUk

· 1

22 Prawdopodobieństwo akceptacji zmiennej jest

równe ilorazowi objętości kuli i opisanej na niej kostki [-1,1]^m

Średnia liczba wylosowanych punktów N_m potrzebnych do realizacji jednej zmiennej wynosi

Modyfikacją usprawniającą powyższy algorytm jest podział kostki na rozłączne podobszary i przeprowadzenia losowania w każdym z nich z osobna.

N

= 1 p

p

= ¼

2

¡(m=2 + 1)

m pm Nm

2 7.854 x 10^-1 1.27 5 1.645 x 10^-1 6.08 10 2.490 x 10^-3 4.015 x 10² 20 2.461 x 10^-8 4.063 x 10⁷ 50 1.537 x 10^-28 6.507 x 10²⁷

Rozkład równomierny na sferze

Jeżeli m-wymiarowa zmienna losowa

ma rozkład równomierny na sferze

oraz macierz A_mxm jest macierzą ortonormalną tj. macierzą obrotu w przestrzeni R^m to zmienna losowa AX ma

również rozkład równomierny na tej sferze.

m-wymiarowa zmienna

ma rozkład sferycznie konturowany jeżeli jej gęstość prawdopodobieństwa

zależy tylko od

X X X = (X

; X

; : : : ; X

)

S

= 8 <

: (x

; : : : ; x

) : X

x

= 1 9 =

;

Z Z Z = (Z

; Z

; : : : ; Z

)

g

(z

; z

; : : : ; z

)

kzk

= X

z

23 Wniosek: mając do dyspozycji odpowiedni rozkład

sferycznie konturowany można go użyć do generowania zmiennej losowej o rozkładzie równomiernym na powierzchnii kuli.

Przykład.

m-wymiarowy rozkład normalny opisuje gęstość

Dla macierzy ortonormalnej zachodzi

więc AZ ma identyczny rozkład.

Algorytm generowania rozkładu równomiernego na sferze w m wymiarach:

1) Generujemy m-wymiarową zmienną losową Z o rozkładzie normalnym

2) Obliczamy

Rozkład nowej zmiennej X będzie równomierny na sferze S_m.

g

(z

; : : : ; z

) = 1

(2¼)

exp µ

¡ 1

2 kzk

¶

kAZ Z Z k

= kZ Z Z k

X

X X = Z Z Z kZ Z Z k

Rozkład równomierny na sferze S_m Tw. Jeżeli zmienna

ma rozkład równomierny na k-wymiarowej sferze S_k, R jest zmienną o rozkładzie

a s jest losowym znakiem

to (k+1)-wymiarowa zmienna losowa

ma rozkład równomierny w (k+1) wymiarowej kuli.

Rozkład równomierny na Kuli K_m

Zauważamy że długość wektora wodzącego R zmiennej X jest także zmienną losową, i ma ona rozkład

Z Z Z = (Z

; Z

; : : : ; Z

)

h(r) =

(

1¡r²

; 0 · r · 1 0 r = 2 [0; 1]

P fs = 1g = P fs = ¡1g = 1 2

X X X = (RZ

; RZ

; : : : ; RZ

; s p

1 ¡ R

)

h(r) = mr

; 0 · r · 1

24 Wystarczy więc wygenerować punkt

leżący na sferze S_m a następnie obliczyć

Zmiena X leży wewnątrz kuli jednostkowej K_m i ma w tym obszarze rozkład równomierny.

Dwuwymiarowy rozkład normalny (metoda Boxa-Mullera)

Jeżeli dwie niezależne zmienne losowe U₁ i U₂ mają rozkład równomierny w przedziale [0,1] to zmienne losowe

również są niezależne a prawdopodobieństwo ich otrzymania opisuje rozkład normalny o

parametrach

Z Z Z = (Z

; Z

; : : : ; Z

)

X

X X = (RZ

; RZ

; : : : ; RZ

)

¹ = 0

¾ = 1 X

= p

¡2ln(U

)cos(2¼U

) X

= p

¡2ln(U

)sin(2¼U

)

Testowanie generatorów liczb pseudolosowych

Ponieważ wszystkie generatory o dowolnym rozkładzie bazują na wykorzystaniu ciągów liczb

losowych o rozkładzie równomiernym więc istotne jest badanie tylko generatorów liczb o takim właśnie

rozkładzie.

Testowanie generatora jest procesem złożonym:

1) Dla ustalonej liczby n, generujemy n kolejny ch liczb startując od losowo wybranej liczby początkowej 2)Obliczamy wartość statystyki testowej (T)

3)Obliczamy F(T) czyli dystrybuantę statystyki T, gdy weryfikowana hipoteza jest prawdziwa

4) Kroki 1-3 powtarzamy N-krotnie obliczając

statystyki: T₁,T₂,...,T_N. Jeśli weryfikowana hipoteza jest prawdziwa to

jest ciągiem zmiennych niezależnych o rozkładzie równomiernym. Testowanie generatora kończy się sprawdzeniem tej hipotezy.

F (T

); F (T

); F (T

); : : : ; F (T

)

25 Testy zgodności z rozkładem równomiernym

U(0,1)

Test chi-kwadrat

Jest najczęściej stosowanym testem.

Badamy w nim hipotezę że generowana zmienna losowa X ma rozkład prawdopodobieństwa o dystrybuancie F.

Jeżeli

to możemy dokonać następującego podziału zbioru wartości zmiennej X

Generujemy n liczb

Sprawdzamy ile z nich spełnia warunek

Ich liczbę oznaczamy n_i.

F (a) = 0 F (b) = 1

a < a

< a

< : : : < a

= b

p

= P fa

< X · a

g; i = 1; 2; : : :

X

; X

; : : : ; X

a

< X · a

Statystyką testu jest

Dla dużego n statystyka ta ma rozkład ² o (k-1) stopniach swobody.

Możemy tak dobrać szerokości przedziałów aby otrzymać zależność

wówczas statystyka przyjmuje prostszą postać

Â

= X

i=1

(n

¡ np

)

np

p

= 1 k

Â

= k n

X

n

¡ n

26 Test OPSO (overlapping-pairs-sparse-occupancy)

Generujemy ciąg liczb

Jeśli z każdej weźmiemy k bitów to możemy utworzyć ciąg liczb całkowitych

z zakresu {0,1,...,2^r-1}.

Następnie tworzymy ciąg kolejnych nakładających się par

Jeśli przez Y oznaczymy liczbę takich par

które nie pojawiły się w ciągu (I_i,I_i+1), to ta zmienna ma rozkład normalny N() – oczywiście dla

odpowiednio dużego n.

Przykładowe parametry testu OPSO

X

; X

; : : : ; X

I

; I

; : : : ; I

(I

; I

); (I

; I

); : : : ; (I

; I

)

f(i; j) : i; j = 0; 1; : : : ; 2

¡ 1g

b n  

10 2

141909 290.26 11 2

1542998 638.75 12 2

567639 580.80

Testy zgodności rozkładów statystyk

Jeżeli wygenerowany ciąg zmiennych losowych

Jest ciągiem zmiennych niezależnych to możemy z elementów tego ciągu utworzyć wektory

które też będą zmiennymi losowymi ale o rozkładzie równomiernym w kostce jednostkowej (0,1)^m.

Możemy zatem zdefiniować pewną funkcję

określoną w kostce jednostkowej.

W ten sposób tworzymy ciąg nowych zmiennych losowych

o jednakowym rozkładzie i dystrybuancie

Testowanie generatora polega na sprawdzeniu hipotezy że ciąg Y₁,Y₂,... jest próbką z populacji o dystrybuancie G.

X

; X

; : : : ; X

(X

; X

; : : : ; X

); (X

; X

; : : : ; X

); : : :

y = h(x

; x

; : : : ; x

)

Y

= h(X

; X

; : : : ; X

) j = 1; 2; : : :

G(y) = P fY

· yg

27 Testy oparte na statystykach

pozycyjnych

Dla wektorów losowych w kostce (0,1)^m definiujemy funkcje

co generuje nowe zmienne

Nowe zmienne mają następujące rozkłady

Testowanie generatora polega na

sprawdzeniu hipotezy o zgodności rozkładów zmiennych U,V,R z powyższym. Testy

przeprowadza się dla niewielkich wartości m=2,3,...,10.

u = max fx

; x

; : : : ; x

g v = min fx

; x

; : : : ; x

g r = u ¡ v

U; V; R

P fU

· ug = u

; 0 · u · 1

P fV

· vg = 1 ¡ (1 ¡ v)

; 0 · v · 1

P fR

· rg = mr

¡ (m ¡ 1)r

; 0 · r · 1

Test sum

Definiujemy funkcję

Wygenerowana zmienna losowa Y ma rozkład o gęstości

Dla m=2

Dla m=3

Testowanie generatora polega na weryfikacji

hipotezy, że zmienna losowa Y ma rozkład zgodny z g_m(y). Zazwyczaj m nie przekracza 5.

y = x

+ x

+ : : : + x

g

(y) =

½ y 0 · y · 1 2 ¡ y 1 < y · 2

g

(y) = 8 <

:

y²

2

0 · y · 1

1

2

(y

¡ 3(y ¡ 1)

) 1 < y · 2

1

2

(y

¡ 3((y ¡ 1)

+ 3(y ¡ 2)

)) 2 < y · 3

28 Testy kombinatoryczne

Pozwalają sprawdzić czy dostarczane przez generator liczby są niezależnymi zmiennymi losowymi.

W testach tych sprawdza się hipotezę że

wielowymiarowa zmienna losowa utworzona z ciągu liczb

ma rozkład o dystrybuancie

Do tej grupy testów zaliczają się np: test pokerowy, test kolekcjonera, test kolizji, test permuatcji.

H(x

; x

; : : : ; x

) = F (x

)F (x

) : : : F (x

)

(X

; X

; : : : ; X

)