Generatory liczb pseudolosowych

(1)

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

(2)

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

0

; X

1

; X

2

; : : : ; X

n

m = 2

³²

U = X

m ) U

ⁱ

2 (0; 1]

X

0

; : : : ; X

º¡1

| {z }

Ta

X

º

; X

º+1

; : : : ; X

º+P¡1

| {z }

To

; X

º+P

; : : :

(3)

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 n.

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

n+1

= (a

1

X

n

+ a

2

X

n¡1

+ : : : +a

k

X

n¡k+1

+ 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

0

; X

1

; X

2

; : : : ; X

k

X

_i+1

= aX

_i_¡1

mod m k

i

=

¹ aX

i¡1

m º

; i ¸ 1 X

1

= aX

0

¡ mk

¹

X

2

= a

²

X

0

¡ mk

²

¡ mk

¹

a

X

3

= a

³

X

0

¡ mk

³

¡ mk

²

a ¡ mk

¹

a

²

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

X

n

= a

ⁿ

X

0

¡ m(k

ⁿ

+ k

n¡1

a + : : : + k

1

a

ⁿ^¡1

)

(4)

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

n

= a

ⁿ

X

0

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

i

= aX

i¡1

mod 11; X

0

= 1

(5)

5 Okres generatora multiplikatywnego

Maksymalny okres generatora multiplikatywnego uzyskujemy dla

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

Przykład

Wykorzystujemy liczby Mersenne'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

0

; i > 0 g

a

^(m^¡1)=p

6= 1 mod m

m = 2

^p

¡ 1

T = 2

³¹

¡ 2

U = fU

ⁱ

; 1 · i · T g

1

2

³¹

¡ 1 = 4:657 £ 10

^¡10

p = 31 ) m = 2

³¹

¡ 1; a = 7

⁵

(6)

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

i

= aX

i¡1

mod 11

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

(U

1

; U

2

; : : : ; U

d

); (U

2

; U

3

; : : : ; U

d+1

); : : :

(U

₁

; U

₂

; : : : ; U

_d

); (U

_d+1

; U

_d+2

; : : : ; U

_2d

); : : :

(7)

7 Parametry statystyczne generatora o rozkładzie

równomiernym w (0,1)-> U(0,1)

Jeśli generowany ciąg 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).

¹ = ¹ 1 N

X

N i=1

x

i

¾

²

= Z

1

0

(x ¡ ¹)

²

dx = 1 12

¾ = ¹ 1 N

X

N i=1

(x

i

¡ ¹ ¹)

²

¹ = Z

1

0

xdx = x

²

2 ¯ ¯

¯ ¯

1

0

= 1

2

(8)

8 Funkcja autokorelacji opisuje zależność

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

Definicja

oraz wzór dla ciągu skończonego (r=s-t)

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.

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

i

= (1176X

i¡1

+ 1476X

i¡2

+ 1776X

i¡3

) mod (2

³²

¡ 5) X

i

= 2

¹³

(X

i¡1

+ X

i¡2

+ X

i¡3

) mod (2

³²

¡ 5)

X

i

= (1995X

i¡1

+ 1998X

i¡2

+ 2001X

i¡3

) mod (2

³⁵

¡ 849) X

i

= 2

¹⁹

(X

i¡1

+ X

i¡2

+ X

i¡3

) mod (2

³²

¡ 1629)

b

i

= (a

1

b

i¡1

+ : : : + a

k

b

i¡k

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

a

1

; a

2

; : : : ; a

k

2 f0; 1g

b

1

; b

2

; : : : ; b

k

2 f0; 1g

R

r

= 1

(N ¡ r)¾

²

N

X

¡r i=1

(X

i

¡ ¹)(X

^i+r

¡ ¹) R ~

r

= E[(X

t

¡ ¹)(X

^s

¡ ¹)]

¾

²

(9)

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

j1

= a

j2

= : : : = a

jk

= 1

b

i

= b

i¡j¹

xor b

i¡j²

xor : : : xor b

i¡j^k

b

i

= b

i¡p

xor b

i¡q

; 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.

Modyfikacja tego typu generatora to generator Mersenne Twister (1998) okres

o rozkładzie

U

i

=

X

L j=1

2

^¡j

b

is+j

= 0:b

is+1

: : : b

is+L

i = 0; 1; 2; : : :

(10)

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

(autokorelacja). Jego modyfikacja

nie posiada już tej wady.

Okres generatora zależy od operacji. Dla

f

0

= f

1

= 1

f

n

= f

n¡2

+ f

n¡1

¦ = +; ¡; ¤; xor

T = 8 >

> <

> >

:

(2

^r

¡ 1)2

^L^¡1

F (r; s; +) (2

^r

¡ 1)2

^L^¡1

F (r; s; ¡) (2

^r

¡ 1)2

^L^¡3

F (r; s; ¤)

(2

^r

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

^L

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

1

; q

2

; : : : ; q

n

)

krk =

Ã

_n

X

i=1

r

_i^p

!

1=p

rrr = (r

1

; r

2

; : : : ; r

n

) P fX = ig = r

ⁱ

P fY = ig = q

ⁱ

i = 1; 2; : : : ; n

±(X) = k(r

¹

; r

2

; : : : ; r

n

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

i

= X

i¡2

+ X

i¡1

mod m; i ¸ 2

X

i

= X

i¡r

¦ X

ⁱ¡s

mod m

i ¸ r; r > s ¸ 1

(11)

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

1

; X

2

; : : : ; X

n

Y

1

; Y

2

; : : : ; Y

n

X

i

¦ Y

^j

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

i

= V

i¡97

¦ V

ⁱ¡33

V

i

2 [0; 1)

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

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

1

; V

2

; : : : ; V

97

V

1

= 0:b

1

b

2

: : : b

24

V

2

= 0:b

25

b

26

: : : b

48

U

n

= V

n

¦ c

ⁿ

(12)

12 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

1

; X

2

; : : : ; X

n

Y

1

; Y

2

; : : : ; Y

n

X

i

¦ Y

^j

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

i

= V

i¡97

¦ V

ⁱ¡33

V

i

2 [0; 1)

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

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

1

; V

2

; : : : ; V

97

V

1

= 0:b

1

b

2

: : : b

24

V

2

= 0:b

25

b

26

: : : b

48

U

n

= V

n

¦ c

ⁿ

(13)

13 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

1

; y

2

; y

3

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

₁

2 f0; 1; : : : ; 168g

c

n

= c

n¡1

¦

µ 7654321 16777216

¶

n ¸ 2

c

1

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

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

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. Wykorzystuje się w nich

„odwrotność modulo”:

Przykład:

Do znalezienia x^-1 wykorzystuje się algorytm podziału Euklidesa (poszukiwania największego wspólnego podzielnika dwóch liczb)

Algorytm:

Jeśli

(gcd-greatest common divisor)

x ¢ x

^¡1

mod m = 1

5 ¢ 3 mod 7 = 1 x = 5 x

^¡1

= 3

a; b; 2 Z; b 6= 0; a = b ¢ q + r q; r 2 Z; 0 · r < jbj

r

_i+1

= 0; = ) gcd(a; b) = r

ⁱ

r

0

= a

r

₁

= b q

_i

=

¹ r

_¡1

r

i

º

.. .

r

i+1

= r

i¡1

¡ q

ⁱ

r

i

; 0 · r

ⁱ⁺¹

< jr

ⁱ

j

(14)

14 Generowanie x^-1 algorytmem Aignera:

1) generujemy ciąg liczb:

i wyznaczamy wspołczynniki q_i

2) generujemy ciąg liczb

3) jeśli x=0 to x^-1=0

Generator nielinowy Eichenauera-Lehna

gdzie: m jest liczbą pierwszą

Uzyskujemy w ten sposób ciąg liczb

o rozkładzie równomiernym.

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

X

i

= (a(i + i

0

) + b)

^¡1

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

z

^m²^¡1

´ 1 mod (z

²

¡ bz ¡ a) fz

ⁿ⁺²

; z

_n+1

; z

_n

; : : : ; z

₁

g

q

_i+1

=

¹ z

i+2

z

_i+1

º

fq

ⁿ⁺¹

; q

_n

; : : : ; q

₁

g

fw

¹

; w

2

; : : : w

n+2

g w

₁

= 1; w

₂

= 0

w

_i

= q

_i_¡1

w

_i_¡1

+ w

_i_¡2

; 3 · i · n + 2 x

^¡1

= ( ¡1)

ⁿ⁺¹

w

_n+2

X

_i+1

= (aX

_i^¡1

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

X

i

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

U

i

= X

i

m 2 [0; 1) z

_n+2

= m; z

_n+1

= x; z

₁

= 1

z

i

= z

i+2

¡ q

ⁱ⁺¹

z

i+1

; 1 · i · n

z 2 Z

(15)

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):

Jeśli uda nam się znaleźć F^-1 to:

zmienna losowa x ma rozkład o dystrybuancie F

F : R ! R

x!¡1

lim F (x) = 0

x

lim

!1

F (x) = 1

U = F (x) ! x = F

^¡1

(U ) F (x) =

Z

x

¡1

f (y)dy

(16)

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

^¡1

(U )

P fX · xg = P fF

^¡1

(U ) · xg

= P fU · F (x)g

= F (x)

U

1

; U

2

; : : : ; U

n

2 (0; 1)

X

1

; X

2

; : : : ; X

n

2 (¡1; 1)

X

1

; X

2

; : : : ; X

n

p

k

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

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

korzystając ze wzoru

Odwracanie dystrybuanty sprawia często duże trudności numeryczne.

U

1

; U

2

; : : : ; U

n

2 (0; 1)

Przykład - rozkład jednomianowy

f (x) = x

ⁿ

n = 1; 2; 3; : : : x 2 [0; 1]

F (y) =

Z

y 0

x

ⁿ

dx = y

ⁿ⁺¹

n + 1 = U

y = ((n + 1)U )

ⁿ⁺¹¹

y 2 [0; 1]

(17)

17 Przykład - rozkład eksponencjalny

gęstość prawdopodobieństwa

Dystrybuanta

f (x) = e

^¡x

; x 2 [0; 1)

F (x) =

Z

x 0

e

^¡x⁰

dx

⁰

F (x) = 1 ¡ e

^¡x

= U

e

^¡x

= 1 ¡ U

U 2 (0; 1) ! x 2 (0; 1)

Przykład - rozkład normalny N(0,1)

F (x) =

Z

x

¡1

e

^¡x⁰²

dx

⁰

= erf (x) f (x) = e

^¡x²

µ

x = y ¡ ¹ p 2¾

¶

(18)

18 Szukanie funkcji odwrotnej erf(x) jest

kosztowne. Częściej stosuje się metodę Boxa-Mullera:

chcemy policzyć prawdopodobieństwo

Wprowadzamy nowe zmienne

f (x; y) = f (x) ¢ f(y) = e

^¡^{x2 +y2}²

p(x; y) = f (x; y)dxdy

r

²

= x

²

+ y

²

x; y 2 (¡1; 1)

r 2 [0; 1) µ 2 [0; 2¼]

p(r; µ) = r ¢ e

^¡r²⁼²

drdµ

Wprowadzamy nową zmienną

z = r

²

2 ! dz = rdr

p(z; µ) = e

^¡z

dzdµ = f (z)dz ¢ dµ

z 2 [0; 1)

f (z) = e

^¡z

Dostajemy rozkład wykładniczy

Ponieważ

µ = U

2

¢ 2¼; U

2

2 (0; 1)

Dla pary (U₁,U₂) dostajemy (x,y) z rozkładu N(0,1)

Przejście

x 2 N(0; 1) ! X 2 N(¹; ¾)

X = x ¢ ¾ + ¹

(19)

19 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.

U

1

2 [a; b] U

2

2 [0; d]

U

2

· f(U

¹

) ) X = U

¹

(20)

20 Generowanie ciągu liczb pseudolosowych o zadanym rozkładzie algorytmem Metropolisa

Modyfikujemy metodę eliminacji. W standardowej postaci jest ona mało wydajna – bo nierzadko odrzucamy większość wyników.

Dzięki algorytmowi Metropolisa nie odrzucamy żadnego.

Akceptację nowego położenia (nowej liczby w ciągu) dokonujemy zgodnie z formułą

czyli:

1) Jeśli f(x_i) > f(x_i-1) to nowe położenie akceptujemy zawsze.

2) W przeciwnym wypadku akceptacja następuje z prawdopodobieństwem f(x_i)/f(x_i-1).

Jeśli nie akceptujemy nowego punktu to zatwierdzamy stary

(każdy krok generuje nowy element ciągu).

Uwaga: w naszym przykładzie

aby zachodził warunek p_ij=p_ji

h(x

_i_¡1

; x

_i

) = min(1; f (x

i

) f (x

_i_¡1

) )

x

_new

= x

_i

+ (2 ¢ U(0; 1) ¡ 1)

x

_i+1

= 8 <

:

x

_i

() x

_new

< 0

x

i

() x

new

> 0 u 2 U(0; 1) > h

x

_new

() x

_new

> 0 u 2 U(0; 1) < h

(21)

21 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ę często zastępuje się sumą

f (x) =

Z

₁

¡1

g

t

(x)h(t)dt

f (x) = X

1

i=1

p

i

g

i

(x) = X

K

i=1

p

i

g

i

(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

i

= Z

A_i

f (x)dx

g

i

(x) = 1 1 1

_A_i

f (x) p

i

I 2 f1; 2; : : : ; Kg

g

i

(x) = 1 1 1

_A_i

f (x) p

i

p

i

¸ 0;

X

K i

p

i

= 1

(22)

22 Przykład – rozkłady o gęstościach

wielomianowych

x 2 [0; 1]

Algorytm generowania zmiennej losowej:

1. Generujemy indeks

z rozkładem prawodopodobieństwa

2. Wylosować zmienną losową X o rozkładzie

Zmienna X ma zadany rozkład wielomianowy.

k 2 f1; 2; 3; : : : ; Mg f (x) =

X

M n=1

c

_n

x

ⁿ

c

n

¸ 0

Z

1 0

f (x)dx = 1 =

X

M n=1

c

_n

n + 1

P fk = ng = c

_n

n + 1

f

_n

(x) = (n + 1)x

ⁿ

Przykład – rozkład Beta Eulera

x 2 [0; 1]

B(p; q) = ¡( p)¡(q)

¡( p + q)

¡( x) =

Z

₁

0

t

^x^¡1

e

^¡t

dt

Jeśli dwie zmienne g₁ i g₂ mają rozkłady

To zmienna Z:

ma rozkład f(x;p,q)

g

1

2 ¡(p) g

2

2 ¡(q)

Z = g

₁

g

₁

+ g

₂

f

_¡

(x; p) = x

^p^¡1

e

^¡x

=¡(p) f (x; p; q) = 1

B(p; q) x

^p^¡1

(1 ¡ x)

^q^¡1

p; q > 0

(23)

23 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

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

X

M k=0

p

k

= 1

U > S

X = X + 1; S = S + p

x

Przykład.

k=0 k=1 k=2 k=3 k=4

U=0.73

0 1.0

p0 = 0.2 p0+p1 = 0.4 p0+p1+p2 = 0.6 p0+p1+p2+p3 = 0.8 p0+p1+p2+p3+p4 = 1

(24)

24 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

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

q

_¡1

= 0 q

i

=

X

i j=0

p

j

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

(25)

25 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

X¡1

> U X = g

X

+ 1

q

X¡1

< U

q

X¡1

> 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_k(0,1)) . Algorytm.

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

Zmienną akceptujemy jeśli

kU U U k

²

· 1 X

X X = (X

1

; X

2

; : : : ; X

k

)

f (x

1

; x

2

; : : : ; x

k

)

U U U = (U

1

; U

2

; : : : ; U

k

)

(26)

26 Prawdopodobieństwo akceptacji zmiennej jest

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

Ś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

m

= 1 p

m

m pk Nk

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 k-wymiarowa zmienna losowa

ma rozkład równomierny na sferze to

k-wymiarowa zmienna X ma rozkład sferycznie konturowany jeżeli jej gęstość prawdopodobieństwa

zależy tylko od

X X X = (X

1

; X

2

; : : : ; X

k

)

S

k

= 8 <

: (x

1

; : : : ; x

k

) : X

k j=1

x

²_j

= 1 9 =

;

g

X

(x

1

; x

2

; : : : ; x

k

) = f ( kxk)

kxk

²

= X

k j=1

x

²_j

(27)

27 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ść

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

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

2) Obliczamy

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

Rozkład równomierny na sferze S_k i S_k+1

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 na (k+1) wymiarowej sferze.

Rozkład równomierny w kuli K_k

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

P fs = 1g = P fs = ¡1g = 1 2 h(r) =

(

pr^k¡1

1¡r²

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

X X X = x x x kxxxk

X

X X = (x

₁

; x

₂

; : : : ; x

_k

)

Z Z Z = (Rx

1

; Rx

2

; : : : ; Rx

k

; s p

1 ¡ R

²

)

(28)

28 Wystarczy więc wygenerować punkt

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

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

Po co rozkład w kuli K_k(0,1)?

Pewne obiekty wielowymiarowe możemy przedstawić w postaci prostych transformacji liniowych współrzędnych.

Taką transformację reprezentuje macierz

Przykład – hiperelipsoidę w R^k otrzymamy transformując K_k(0,1).

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

1

); F (T

2

); F (T

3

); : : : ; F (T

N

) Z

Z Z = (Z

₁

; Z

₂

; : : : ; Z

_k

)

X X X = (RZ

1

; RZ

2

; : : : ; RZ

k

)

Z = A ~ ~ X

(29)

29 Testy zgodności z zadanym rozkładem

Test chi-kwadrat

(rozkład ciągły lub dyskretny) 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

1

< a

2

< : : : < a

k

= b

p

i

= P fa

ⁱ¡1

< X · a

ⁱ

g; i = 1; 2; : : :

X

1

; X

2

; : : : ; X

n

a

i¡1

< 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ć

Â

²_k_¡1

= X

k

i=1

(n

i

¡ np

ⁱ

)

²

np

i

p

i

= 1 k

Â

²_k_¡1

= k n

X

k i=1

n

²_i

¡ n

(30)

30 Przykład. Test generatora o trójkątnym rozkładzie fgp.

Test -procedura:

1. Generujemy serię N liczb z rozkładu 2. Grupujmey liczby w k przedziałach

3. Liczymy wartość statystyki testowejchi^2 4. Dla zadanego poziomu  sprawdzamy czy

chi^2 < wartości granicznej

(31)

31 Test możemy rozszerzyć – test zgodności z rozkładem ciągłym Kołmogorova-Smirnova

Generujemy M próbek o rozkładzie trójkątnym (N liczb w każdej próbce) Dla każdej próbki obliczamy wartość statystyki testowej ² i sortujemy je.

Dla każdej wartości obliczamy empiryczną wartość dystrybuanty i jej dokładną wartość lub grupujemy je tak aby częstości wystęþowania wzrastające o np.: p=0.05

Obliczamy wartości statystyk

Ich wartości porównujemy z wartościami granicznymi dla założonego poziomu ufności, odrzucając lub nie hipotezę że wygenerowany ciąg liczb ma zadany (np. trójkątny) rozkład.

(32)

32 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

1

; X

2

; : : : ; X

n

I

1

; I

2

; : : : ; I

n

(I

1

; I

2

); (I

2

; I

3

); : : : ; (I

n¡1

; I

n

)

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

^b

¡ 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

1

; X

2

; : : : ; X

n

(X

1

; X

2

; : : : ; X

m

); (X

m+1

; X

m+2

; : : : ; X

2m

); : : :

y = h(x

1

; x

2

; : : : ; x

m

)

Y

j

= h(X

_(j_¡1)m+1

; X

_(j_¡1)m+2

; : : : ; X

jm

) j = 1; 2; : : :

G(y) = P fY

^j

· yg

(33)

33 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.

r = u ¡ v

U; V; R

P fU

^j

· ug = u

^m

; 0 · u · 1

P fV

^j

· vg = 1 ¡ (1 ¡ v)

^m

; 0 · v · 1

P fR

^j

· rg = mr

^m^¡1

¡ (m ¡ 1)r

^m

; 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

1

+ x

2

+ : : : + x

m

g

2

(y) =

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

g

3

(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