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 RM 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 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
Ta – okres aperiodyczności ciągu T0 – 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
nm = 2
32U = X
m ) U
i2 (0; 1]
X
0; : : : ; X
º¡1| {z }
Ta
X
º; X
º+1; : : : ; X
º+P¡1| {z }
To
; X
º+P; : : :
3 Generatory liniowe
Generatory liniowe tworzą ciąg liczb według schematu:
gdzie:
a1, a2,...,ak, 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
1X
n+ a
2X
n¡1+ : : : +a
kX
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 (X0).
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
kX
i+1= aX
i¡1mod m k
i=
¹ aX
i¡1m
º
; i ¸ 1 X
1= aX
0¡ mk
1X
2= a
2X
0¡ mk
2¡ mk
1a
X
3= a
3X
0¡ mk
3¡ mk
2a ¡ mk
1a
2: : : : : : : : : : : :
X
n= a
nX
0¡ m(k
n+ k
n¡1a + : : : + k
1a
n¡1)
4 Ostatnie równanie można zapisać w postaci
skąd wynika, że wybór X0 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
nX
0mod 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¡1mod 11; X
0= 1
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 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
i= X
0; i > 0 g
a
(m¡1)=p6= 1 mod m
m = 2
p¡ 1
T = 2
31¡ 2
U = fU
i; 1 · i · T g
1
2
31¡ 1 = 4:657 £ 10
¡10p = 31 ) m = 2
31¡ 1; a = 7
56 Rozkład przestrzenny ciągu
Wadą generatorów multiplikatywnych jest
nierównomierne pokrycie d-wymiarowej kostki (Id).
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¡1mod 11
a) X0=1, a=2 b) X0=1, a=8
(U
1; U
2; : : : ; U
d); (U
2; U
3; : : : ; U
d+1); : : :
(U
1; U
2; : : : ; U
d); (U
d+1; U
d+2; : : : ; U
2d); : : :
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=1x
i¾
2= Z
10
(x ¡ ¹)
2dx = 1 12
¹
¾ = 1 N
X
N i=1(x
i¡ ¹¹)
2Przykł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
i+1= aX
imod 2
Lc ´ 3 mod 8; X
0´ 1; 3; 9; 11 mod 16 X
0; X
1; X
2; : : :
8j + 1; 8j + 3; j = 0; 1; : : : ; 2
L¡3¡ 1 T = 2
L¡2¾
2= 1
12 ¡ 13 3 ¢ 2
m¹
¹ = 1
2 ¡ 1 2
L¡1¹ = Z
10
xdx = x
22
¯ ¯
¯ ¯
1 0
= 1 2
U
i= X
im = X
i2
L8 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.
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
12 [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 bi otrzymywanych rekurencyjnie
gdzie:
są stałymi binarnymi , a stałe bi
tworzą ciąg inicjujący.
X
i= (1176X
i¡1+ 1476X
i¡2+ 1776X
i¡3) mod (2
32¡ 5) X
i= 2
13(X
i¡1+ X
i¡2+ X
i¡3) mod (2
32¡ 5)
X
i= (1995X
i¡1+ 1998X
i¡2+ 2001X
i¡3) mod (2
35¡ 849) X
i= 2
19(X
i¡1+ X
i¡2+ X
i¡3) mod (2
32¡ 1629)
b
i= (a
1b
i¡1+ : : : + a
kb
i¡k) mod 2 i = k + 1; k + 2; : : :
a
1; a
2; : : : ; a
k2 f0; 1g
b
1; b
2; : : : ; b
k2 f0; 1g
R
r= 1
(N ¡ r)¾
2N
X
¡r i=1(X
i¡ ¹)(X
i+r¡ ¹) R ~
r= E[(X
t¡ ¹)(X
s¡ ¹)]
¾
29 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 2k a dokładniej (2k-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¡j1xor b
i¡j2xor : : : xor b
i¡jkb
i= b
i¡pxor 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 Ui oraz Ui+1
wykorzystywane są elementy tego samego podciągu
b) jeśli s=L to Ui oraz Ui+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
i=
X
L j=12
¡jb
is+j= 0:b
is+1: : : b
is+Li = 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
(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¡1F (r; s; +) (2
r¡ 1)2
L¡1F (r; s; ¡) (2
r¡ 1)2
L¡3F (r; s; ¤)
(2
r¡ 1) F (r; s; xor) m = 2
LGenerator 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 =
Ã
nX
i=1
r
ip!
1=prrr = (r
1; r
2; : : : ; r
n) P fX = ig = r
iP fY = ig = q
ii = 1; 2; : : : ; n
±(X) = k(r
1; r
2; : : : ; r
n) ¡ (1=n; 1=n; : : : ; 1=n)k X
i= X
i¡2+ X
i¡1mod m; i ¸ 2
X
i= X
i¡r¦ X
i¡smod 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 T1, a ciąg
ma okres T2 i są liczbami pierwszymi to wtedy nowy ciąg liczb
ma okres równy T1T2.
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
nY
1; Y
2; : : : ; Y
nX
i¦ Y
jPierwszy 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
i¡33V
i2 [0; 1)
x ¦ y = x ¡ y; x ¸ y F (97; 33; ¦)
x ¦ y = x ¡ y + 1; x < y V
1; V
2; : : : ; V
97V
1= 0:b
1b
2: : : b
24V
2= 0:b
25b
26: : : b
48y
n= (y
n¡3y
n¡2y
n¡1) mod 179 z
n= (52z
n¡1+ 1) mod 169 b
n=
½ 0; gdy y
nz
nmod 64 < 32
b
n= 1; w pozostalych przypadkach
U
n= V
n¦ c
n12 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 T1, a ciąg
ma okres T2 i są liczbami pierwszymi to wtedy nowy ciąg liczb
ma okres równy T1T2.
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
nY
1; Y
2; : : : ; Y
nX
i¦ Y
jPierwszy 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
i¡33V
i2 [0; 1)
x ¦ y = x ¡ y; x ¸ y F (97; 33; ¦)
x ¦ y = x ¡ y + 1; x < y V
1; V
2; : : : ; V
97V
1= 0:b
1b
2: : : b
24V
2= 0:b
25b
26: : : b
48y
n= (y
n¡3y
n¡2y
n¡1) mod 179 z
n= (52z
n¡1+ 1) mod 169 b
n=
½ 0; gdy y
nz
nmod 64 < 32
b
n= 1; w pozostalych przypadkach
U
n= V
n¦ c
n13 Generatory muszą być zainicjowane
Okres generatora kombinowanego: 2120.
Drugi generator (cn) o rozkładzie równomiernym w (0,1) jest zdefiniowany następująco
Okres tego generatora to 224-3.
Okres generatora uniwersalnego (Un): 2144
y
1; y
2; y
32 f1; 2; : : : ; 178g z
12 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
¡1mod 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
ir
0= a
r
1= b q
i=
¹ r
¡1r
iº
.. .
r
i+1= r
i¡1¡ q
ir
i; 0 · r
i+1< jr
ij
14 Generowanie x-1 algorytmem Aignera:
1) generujemy ciąg liczb:
i wyznaczamy wspołczynniki qi
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
2¡ 1
X
i= (a(i + i
0) + b)
¡1mod m; i = 0; 1; : : :
z
m2¡1´ 1 mod (z
2¡ bz ¡ a) fzn+2; z
n+1; z
n; : : : ; z
1g
q
i+1=
¹ z
i+2z
i+1º
fq
n+1; q
n; : : : ; q
1g fw
1; w
2; : : : w
n+2g w
1= 1; w
2= 0
w
i= q
i¡1w
i¡1+ w
i¡2; 3 · i · n + 2 x
¡1= ( ¡1)
n+1w
n+2X
i+1= (aX
i¡1+ b) mod m; i = 0; 1; : : :
X
i2 f0; 1; : : : ; m ¡ 1g
U
i= X
im 2 [0; 1) zn+2 = m; z
n+1 = x; z
1 = 1
z
i= z
i+2¡ q
i+1z
i+1; 1 · i · n
z 2 Z
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
!1F (x) = 1
U = F (x) ! x = F
¡1(U ) F (x) =
Z
x¡1
f (y)dy
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 Xi 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
n2 (0; 1)
X
1; X
2; : : : ; X
n2 (¡1; 1)
X
1; X
2; : : : ; X
np
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
n2 (0; 1)
X
n= min (
k : U
n· X
ki=0
p
i)
; n = 1; 2; : : :
Przykład - rozkład jednomianowy
f (x) = x
nn = 1; 2; 3; : : : x 2 [0; 1]
F (y) = Z
y0
x
ndx = y
n+1n + 1 = U
y = ((n + 1)U )
n+11y 2 [0; 1]
17 Przykład - rozkład wykładniczy
gęstość prawdopodobieństwa
Dystrybuanta
f (x) = e
¡x; x 2 [0; 1)
F (x) =
Z
x 0e
¡x0dx
0F (x) = 1 ¡ e
¡x= U e
¡x= 1 ¡ U
F
¡1(x) = x = ¡ln(1 ¡ U) U 2 (0; 1) ! x 2 (0; 1)
Przykład - rozkład normalny N(0,1)
F (x) =
Z
x¡1
e
¡x02dx
0= erf (x) f (x) = e
¡x2µ
x = y ¡ ¹ p 2¾
¶
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 +y22p(x; y) = f (x; y)dxdy
r
2= x
2+ y
2x = rcos(µ)
y = rsin(µ)
x; y 2 (¡1; 1)
r 2 [0; 1) µ 2 [0; 2¼]
p = f (x; y)dxdy = f (r; µ)drdµ p(r; µ) = r ¢ e
¡r2=2drdµ
Wprowadzamy nową zmienną
z = r
22 ! dz = rdr
p(z; µ) = e
¡zdzdµ = f (z)dz ¢ dµ
z 2 [0; 1)
f (z) = e
¡zDostajemy rozkład wykładniczy
z = ¡ln(1 ¡ U
1); U
12 (0; 1)
x = rcos(µ) = p
¡2ln(1 ¡ U
1)cos(2¼U
2) y = rsin(µ) = p
¡2ln(1 ¡ U
1)sin(2¼U
2)
Ponieważ
µ = U
2¢ 2¼; U
22 (0; 1)
Dla pary (U1,U2) dostajemy (x,y) z rozkładu N(0,1)
Przejście
x 2 N(0; 1) ! X 2 N(¹; ¾)
X = x ¢ ¾ + ¹
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ę U1, U2
4) wykonujemy czynności 1-3 aż do uzyskania odpowiednio licznego ciągu
Wygenerowana zmienna losowa X ma rozkład prawdopodobieństwa f.
U
12 [a; b] U
22 [0; d]
U
2· f(U
1) ) X = U
120 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(xi) > f(xi-1) to nowe położenie akceptujemy zawsze.
2) W przeciwnym wypadku akceptacja następuje z prawdopodobieństwem f(xi)/f(xi-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 pij=pji
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 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: gt(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 gt(x) W praktyce całkę często zastępuje się sumą
f (x) =
Z
1¡1
g
t(x)h(t)dt
f (x) = X
1i=1
p
ig
i(x) = X
Ki=1
p
ig
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 (Ai) wyznaczamy
oraz
gdzie: 1Ai jest funkcją przynależności do podzbioru Ai. 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
Ai
f (x)dx
g
i(x) = 111
Aif (x) p
iI 2 f1; 2; : : : ; Kg
g
i(x) = 111
Aif (x) p
ip
i¸ 0;
X
K ip
i= 1
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=1c
nx
nc
n¸ 0 Z
10
f (x)dx = 1 =
X
M n=1c
nn + 1
P fk = ng = c
nn + 1
f
n(x) = (n + 1)x
nPrzykład – rozkład Beta Eulera
x 2 [0; 1]
B(p; q) = ¡(p)¡(q)
¡(p + q)
¡(x) =
Z
10
t
x¡1e
¡tdt
Jeśli dwie zmienne g1 i g2 mają rozkłady
To zmienna Z:
ma rozkład f(x;p,q)
g
12 ¡(p) g
22 ¡(q)
Z = g
1g
1+ g
2f
¡(x; p) = x
p¡1e
¡x=¡(p) f (x; p; q) = 1
B(p; q) x
p¡1(1 ¡ x)
q¡1p; q > 0
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=p0
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 X1,X2,X3,...
P fX = kg = p
k; k = 0; 1; 2; : : : ; M
X
M k=0p
k= 1
U > S
X = X + 1; S = S + p
xb) 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
g
i= max
½
j : q
j< i K + 1
¾
; i = 1; 2; : : : ; K + 1 q
i=
X
i j=0p
j; i = 0; 1; : : : ; K
24 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 X1,X2,...
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 (Kk(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
2· 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)
25 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 Nm 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
mm 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 102 20 2.461 x 10-8 4.063 x 107 50 1.537 x 10-28 6.507 x 1027
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=1x
2j= 1 9 =
;
g
X(x
1; x
2; : : : ; x
k) = f ( kxk)
kxk
2= X
k j=1x
2j26 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 Sk.
Rozkład równomierny na sferze Sk i Sk+1
Tw. Jeżeli zmienna
ma rozkład równomierny na k-wymiarowej sferze Sk, 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 Kk
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) = (2k ¡ 1)r
k¡1; 0 · r · 1 h(r) =
(
prk¡11¡r2
; 0 · r · 1 0 r = 2 [0; 1]
X X X = x x x kxxxk
X X X = (x
1; x
2; : : : ; x
k)
Z Z Z = (Rx
1; Rx
2; : : : ; Rx
k; s p
1 ¡ R
2) g
X(x
1; : : : ; x
k) = 1
(2¼)
k=2exp µ
¡ 1
2 k~xk
2¶
27 Wystarczy więc wygenerować punkt
leżący na sferze Sk a następnie obliczyć
Zmiena X leży wewnątrz kuli jednostkowej Kk i ma w tym obszarze rozkład równomierny.
Po co rozkład w kuli Kk(0,1)?
Pewne obiekty wielowymiarowe możemy przedstawić w postaci prostych transformacji liniowych współrzędnych.
Taką transformację reprezentuje macierz
Przykład – hiperelipsoidę w Rk otrzymamy transformując Kk(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: T1,T2,...,TN. 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
1; Z
2; : : : ; Z
k)
X X X = (RZ
1; RZ
2; : : : ; RZ
k)
Z = A ~ ~ X
28 Testy zgodności z zadanym rozkładem
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 ni.
F (a) = 0 F (b) = 1
a < a
1< a
2< : : : < a
k= b
p
i= P fa
i¡1< X · a
ig; i = 1; 2; : : :
X
1; X
2; : : : ; X
na
i¡1< X · a
iStatystyką testu jest
Dla dużego n statystyka ta ma rozkład 2 o (k-1) stopniach swobody.
Możemy tak dobrać szerokości przedziałów aby otrzymać zależność
wówczas statystyka przyjmuje prostszą postać
Â
2k¡1= X
ki=1
(n
i¡ np
i)
2np
ip
i= 1 k
Â
2k¡1= k n
X
k i=1n
2i¡ n
29 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,...,2r-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 (Ii,Ii+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
nI
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
21141909 290.26 11 2
221542998 638.75 12 2
23567639 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 Y1,Y2,... 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
30 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
1; x
2; : : : ; x
mg v = min fx
1; x
2; : : : ; x
mg 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 gm(y). Zazwyczaj m nie przekracza 5.
y = x
1+ x
2+ : : : + x
mg
2(y) =
½ y 0 · y · 1 2 ¡ y 1 < y · 2
g
3(y) = 8 <
:
y2
2
0 · y · 1
1
2
(y
2¡ 3(y ¡ 1)
2) 1 < y · 2
1
2