Komputerowe generatory liczb losowych [PDF]

Zitiervorschau

VN1 \ V

7

\ *\ \

Spis treści Przedmowa...................................................................................................................................................... Wykaz niektórych oznaczeń............................................................................................................................. 1. „ Liczby losowe".......................................................................................................................................... 2. Generatory liczb losowych o rozkładzie równomiernym.............................................................................. 2.1. Wprowadzenie....................................................................................................................................... 2.2. Generatory liniowe................................................................................................................................ 2.2.1. Opis.................................................................................................................................................... 2.2.2. Okres generatora............................................................................................................................. 2.2.3. Struktura przestrzenna.................................................................................................................... 2.2.4. Ogólne generatory liniowe.............................................................................................................. 2.2.5. Parametry statystyczne................................................................................................................... 2.2.6. Wybór parametrów dla generatorów liniowych............................................................................... 2.3. Generatory oparte na rejestrach przesuwnych........................................................................................ 2.4. Generatory Fibonacciego....................................................................................................................... 2.5. Kombinacje generatorów....................................................................................................................... 2.6. Uniwersalny generator liczb losowych o rozkładzie równomiernym...................................................... 2.7. Generatory oparte na odejmowaniu z pożyczką i generator ULTRA...................................................... 2.8. Generatory oparte na mnożeniu z przeniesieniem.................................................................................. 2.9. Generatory nieliniowe............................................................................................................................ 2.10. Uwagi o implementacji numerycznej................................................................................................... 2.11. Przykładowe implementacje w języku C .............................................................................................. 2.11.1. Inicj owanie generatorów....................................................................................................... 2.11.2. Ogólny generator liniowy..................................................................................................... 2.11.3. Generator Tauswortha z podrozdziału 2.3..................................................................................... 2.11.4. Generator uniwersalny z podrozdziału 2.6..................................................................................... 3. Generatory liczb losowych o dowolnych rozkładach prawdopodobieństwa.................................................. 3.1. Ogólne metody konstrukcji generatorów liczb losowych o dowolnych rozkładach prawdopodobieństwa 3.1.1. Metoda odwracania dystrybuanty.................................................................................................... 3.1.2. Metoda eliminacji........................................................................................................................... 3.1.3. Metoda superpozycji rozkładów...................................................................................................... 3.1.4. Metoda ROU................................................................................................................................... 3.1.5. Rozkłady dyskretne......................................................................................................................... 3.2. Metody konstrukcji generatorów dla podstawowych rozkładów prawdopodobieństwa.................... 3.2.1. Rozkłady dyskretne......................................................................................................................... 3.2.2. Rozkład wykładniczy...................................................................................................................... 3.2.3. Rozkład normalny........................................................................................................................... 3.2.4. Rozkład gamma.............................................................................................................................. 3.2.5. Rozkład beta................................................................................................................................... 3.2.6. Rozkład Cauchy'ego....................................................................................................................... 3.2.7. Rozkłady a-stabilne........................................................................................................................ 3.3. Związki między rozkładami................................................................................................................... 4. Generatory liczb losowych o rozkładach wielowymiarowych....................................................................... 4.1. Przypadek ogólny.................................................................................................................................. 4.2. Rozkłady równomierne w Rm................................................................................................................. 4.2.1. Uwagi ogólne.................................................................................................................................. 4.2.2. Rozkład równomierny na sferze i na kuli w Rm............................................................................... 4.2.3. Rozkład równomierny na sympleksie i na powierzchni sympleksu.................................................. 4.3. Wielowymiarowy rozkład normalny...................................................................................................... 5. Testowanie generatorów liczb losowych....................................................................................................... 5.1. Metodyka testowania generatorów......................................................................................................... 5.2. Testy zgodności z rozkładem U(0,1)...................................................................................................... 5.2.1. Test chi-kwadrat............................................................................................................................. 5.2.2. Test zgodności z rozkładem wielowymiarowym............................................................................. 5.2.3. Test OPSO...................................................................................................................................... 2.4. Test Kołmogorowa................................................................................................................................ 5.3. Testy zgodności rozkładów statystyk..................................................................................................... 5.3.1. Wprowadzenie................................................................................................................................ 5.3.2. Testy oparte na statystykach pozycyjnych....................................................................................... 5.3.3. Test sum.......................................................................................................................................... 5.3.4. Test d2............................................................................................................................................ 5.3.5. Test urodzin dla spacji....................................................................................................................

.. 4 .. 5 .. 6 .. 8 .. 8 10 10 11 12 13 13 14 16 19 20 21 22 23 24 26 26 26 27 27 28 30 30 30 32 43 50 53 56 56 61 63 67 72 74 76 77 82 82 83 83 84 86 88 89 89 91 91 91 92 93 94 95 95 96 96 97

5.3.6. Test najmniejszej odległości w parach...................................................................................................... 97 5.5. Testy serii....................................................................................................................................................... 100 5.6. Testy kombinatoryczne.................................................................................................................................... 102 5.6.1. Testpokerowy..........................................................................................................................................102 5.6.2. Testkolekcjonera...................................................................................................................................... 104 5.6.3. Test kolizji i test liczby pustych cel...........................................................................................................104 5.6.4. Testpermutacji.........................................................................................................................................105 5.6.5. Test oparty na rzędzie losowych macierzy binarnych................................................................................105 5.7. Testowanie generatorów za pomocą zadań kontrolnych................................................................................. 105 6.Prace cytowane....................................................................................................................................................... 106

Przedmowa Losow anie prób w kontekście badań statystycznych (badania reprezentacyjne, sym ulacyjne badania estym atorów , testów i statystycznych reguł decyzyjnych) oraz w kontekście obliczeń numerycznych (metody M onte Carlo, klasyczne dla całek i rów nań z operatorami liniowymi i now sze dla zadań optym alizacji), ja k rów nież sym ulacyjne badania m odeli probabilistycznych w technice, ekonomii, naukach przyrodniczych i praktycznie we w szystkich dziedzinach wiedzy, wymagają wyposażenia współczesnego kom putera w odpow iednie narzędzia. Takimi narzędziam i są generatory liczb losowych. K om putery trafiły pod strzechy i fala publikacji poświęconych różnym aspektom ich w y k o rzy stan ia nie opada. Jesteśm y p rzek o n an i, że w czasie, jaki upłynie m iędzy postaw ieniem przez nas ostatniej kropki w kom puteropisie tej książki a jej dotarciem do pierw szych Czytelników , pojawi się co najmniej kilkadziesiąt nowych publikacji i program ów bezpośrednio zw iązanych z tem atyką generatorów liczb losow ych. W ierzym y jednak, że to co proponujem y Czytelnikom, będzie jeszcze przez pew ien czas aktualne, a przynajm niej ułatw i Im, i jeszcze długo będzie ułatw iało, poruszanie się w gąszczu coraz to now ych osiągnięć w tej dziedzinie. W arszaw a, maj 1997

A utorzy

Wykaz niektórych oznaczeń

T

a [a] {a} Ml a mod b a xor b A C ov(X , Y) D2(X) E (X ) 1a lm(A) In a P{A} R Rn sign m o rn Q □

transpozycja w ektora a n a jw ię k sz a lic z b a całk o w ita, m n ie jsz a lub ró w n a a ułam kow a część liczby a norm a (długość) w ektora a reszta z dzielenia liczby a p rze z b operacja binarna (a + b) m od 2 dopełnienie zbioru A kow ariancja zm iennych losow ych X i Y w ariancja zm iennej losowej X w artość oczekiw ana zm iennej losowej X funkcja charakterystyczna zbioru A m iara (L ebesgue'a) zbioru A w R m logarytm naturalny liczby a praw dopodobieństw o zdarzenia A zbiór liczb rzeczywistych przestrzeń euklidesow a n-w ym iarow a funkcj a znaku funkcja gam m a E ulera dystrybuanta rozkładu norm alnego A (0,1) zdarzenie elem entarne przestrzeń zdarzeń elem entarnych koniec przykładu koniec lem atu lub tw ierdzenia

Any one who considers arithmetical methods ofpro-ducing random digits is, ofcourse, in a state ofsin. For, as has been pointed out severa/ times, there is no suci thing as a random number - there are only methods to pro-duce random numbers, and a strict arithmetic procedure ofcourse is not such

a method.1'

JOHN VON NEUMANN, 1951

1. „ Liczby losowe” W ykonajm y serię niezależnych rzutów m onetą i zanotujmy obserwowane w yniki, pisząc, O gdy w ynikiem rzutu je st reszka, lub l - gdy w ynikiem rzutu je st orzeł. Poniew aż praw dopodobieństw o zaobserw ow ania orla je st takie sam o ja k praw dopodobieństw o zaobserw ow ania reszki i rów ne 1/2, w ięc zm ienna losow a w ynik rzutu m a rozkład dwupunktow y i przyjm uje wartości O lub l z jednakowym prawdopodobieństwem. Mówimy, że ta zm ienna losowa m a rozkład równomierny na zbiorze {0,1}. W wyniku opisanego postępow ania otrzym am y, w ięc np. następujący ciąg liczb: 1 ,0 ,0 ,1 ,1 ,... Taki ciąg przyjęto nazywać ciągiem liczb losowych o rozkładzie równomiernym na zbiorze {0,1}. M onetę, za pom ocą, której otrzym ujem y tego typu ciągi, nazywamy generatorem liczb losowych o rozkładzie równomiernym na zbiorze {0,1}. P rzygotujm y 10000 jed n ak o w y ch kartek i ponum erujm y je kolejnym i liczbami czterocyfrowymi: 0000,0001,0002,..., 9999. W rzućmy wszystkie kartki do urny. Losujm y z urny po jednej kartce, tak żeby w każdym losow aniu k ażda z nich m iała jed n ak o w ą szansę na w yjęcie z urny. Z m ienna losow a w ylosow any num er m a ro zkła d rów nom ierny na zbiorze {0000,0001,0002,

1) Jawnie grzeszy, kto opowiada o arytmetycznych procedurach generowania liczb losowych. Nie ma bowiem, jak to już nieraz mówiono, czegoś takiego, jak liczba losowa. Istnieją metody losowego wytwarzania liczb, ale oczywiście żadna deterministyczna procedura arytmetyczna nie jest taką metodą. (Tłumaczenie autora)

9999}. Jeżeli opisane losowanie pow tórzym y w ielokrotnie (po każdymi z nich zw racając w ylosow aną kartkę z powrotem do urny), to otrzym amy np. następujący ciąg liczb: 1722,4355, 0234,... Taki ciąg będziem y nazywać t ciągiem liczb losowych o rozkładzie równomiernym na zbiorze {0000,0001,1 0002,..., 9999}, a urnę z ponumerowanym i kartkami - generatorem liczb losowych o rozkładzie równomiernym na tym zbiorze. W eźm y pod uw agę ruletkę z kołem o obwodzie równym l i niech A będzie w yróżnionym punktem na obwodzie tarczy tej ruletki. Tarczę wprawiaj się w ruch obrotowy i rzuca na nią kulkę K, która, dzięki odpowiedniej konstrukcji ruletki, zatrzym uje się zaw sze na brzegu tarczy. N iech U będzie długością łuku AK. Załóżmy, że zm ienna losow a U m a rozkład rów nom ierny na przedziale (0,1) - rozkład ten będziem y oznaczali przez U(0,1). W wyniki kilku eksperym entów z ruletką otrzym amy np. następujący ciąg liczb: 0.2217,1 0.5543, 0.3402,... Taki ciąg przyjęto nazywać ciągiem liczb losowych o rozkładzie równomiernym na przedziale (0,1), a ruletkę - generatorem liczb losowych o rozkładzie równomiernym U(0,1). Zadania, w których do rozw iązywania używ a się ciągów liczb losowych, m ożna podzielić na trzy grupy. Grupę pierw szą (również pierw szą historycznie) tw orzą zadania związanej z badaniami reprezentacyjnymi. Problem opisu różnych zbiorów (populacji) za pom ocą próbek losow anych z tych zbiorów je st typow ym problem em statystycznym. Przykładam i są tu badania różnych zjawisk społecznych przez szczegółowy opis jednostek w ybranych losow o z populacji interesujących badacza obiektów (ludzi, zakładów pracy, środowisk, szkół, itp.) lub zadania ze statystycznej kontroli jakości, w których partie różnych tow arów opisuje się na podstaw ie badania losowo w ybranych próbek tych towarów. W praktyce stosuje się tu nie tylko takie schematy losowania, jak podany wyżej przykład losow ania prostego z urny; obszerny przegląd różnych m etod generow ania prób losow ych w takich sytuacjach m ożna znaleźć np. w książce Zasępy (1962) i w książce Brachy (1996). Grupę drugą stanowią zadania num eryczne rozw iązywane metodami M onte Carlo. Zadanie num eryczne (typowym przykładem je st zadanie obliczania wartości danej całki) zastępuje się w ówczas zadaniem rachunku prawdopodobieństwa, które z kolei rozw iązuje się na drodze eksperym entu statystycznego. Podstaw ow ą częścią tego eksperym entu je st losowanie próbki z odpowiedniej populacji, a w ięc generow anie odpowiedniego ciągu liczb losowych. Obszerny w ykład różnych sposobów postępow ania w takich sytuacjach m ożna znaleźć w książkach: H am m ersleya i H andscom ba (1961), Zielińskiego (1970), Jerm akow a (1976), N iederreitera i Shiue (1995). N ajnow sze metody tego typu dotyczą stochastycznych algorytm ów szukania minimum globalnego danej funkcji, czemu pośw ięcona je st m inim onografia Zielińskiego i N eum anna (1986) oraz liczne prace dotyczące sym ulowanego w yżarzania (ang. sim ulated annealing) z ostatniego dziesięciolecia; aktualne wyniki na ten tem at m ożna znaleźć w pracy W ieczorkow skiego (1995). Grupę trzecią stanowią zadania zw iązane z badaniem różnych zjaw isk i procesów (technicznych, ekonom icznych, przyrodniczych) za pom ocą ich komputerowej symulacji (modelowania). O przebiegu takich procesów de-dują najczęściej czynniki losowe, a m odelowanie w pływ u tych czynników sprowadza się do losow ania próbek z odpowiednich rozkładów prawdopodobieństwa, czyli do generow ania odpow iednich ciągów liczb losowych. W sytuacjach realnych korzystanie z opisanych wyżej generatorów liczb losowych (moneta, urna lub ruletka) jest, oczywiście, najczęściej niem ożliwe j w praktyce takie „prawdziwe" generatory są zastępowane pewnymi ich nam iastkam i. Jeszcze do niedaw na powszechnie posługiw ano się tablicami liczb losowych, obecnie stosuje się odpowiednie programy komputerowe. W szędzie dalej w tej książce m ówiąc o generatorach liczb losowych mamy właśnie na myśli takie programy. Podstaw ow ą rolę odgryw ają generatory liczb losowych o rozkładzie równomiernym U (0,1) - omawiamy je w rozdz. 2. Liczby otrzym ywane w w yniku obliczeń wykonywanych za pom ocą program ów kom puterowych nie są oczywiście „tak losowe" ja k liczby uzyskiw ane przez rzuty monetą, losowanie z urny lub obracanie koła ruletki. W celu podkreślenia tego faktu używ a się często nazw liczby pseudolosow e lub liczby ąuasi-losowe, ale nie będziem y

tutaj rygorystycznie trzym ali generatorach liczb losowych.

się tych nazw;

dalej

m ów im y po prostu o program ow ych

Jak ju ż wspom nieliśmy, w zastosow aniach potrzebne są liczby losow e o różnych rozkładach prawdopodobieństwa, np. liczby o rozkładzie norm alnym lub liczby opisujące realizacje procesu Poissona. W szystkie takie liczby losow e m ożna otrzym ać przez odpow iednie m anipulacje liczbami z generatora liczb losowych o rozkładzie równomiernym U(0,1). M ówim y o tym dokładnie w rozdz. 3. (przypadek rozkładów jednow ym iarow ych) i w rozdz. 4. (rozkłady wielowym iarowe). W zw iązku ze stosowaniem różnych generatorów liczb losow ych pow staje problem testowania tych generatorów. Ogólnie mówiąc, sprowadza się on do testow ania odpowiednich hipotez statystycznych o generatorze, co omówimy w rozdz. 5. Kom puterow e generatory liczb losowych są w dzisiejszych czasach jednym z najczęściej używ anych narzędzi każdego badacza: matematyka, lekarza, inżyniera, ekonomisty, socjologa, chem ika i fizyka - do symulacji procesów losowych. U żytkow nik tego kom puterowego narzędzia zwykle bardzo szybko zaczyna doceniać jego potężne m ożliwości, a zarazem, posługując się nim, odczuw a przyjem ność i satysfakcję. Życzym y tego naszym Czytelnikom.

2. Generatory liczb losowych o rozkładzie równomiernym

2.1. Wprowadzenie N ajprostszym i generatorami liczb losowych są oczywiście generatory fizyczne, ja k np. w ym ienione w rozdz. 1. moneta, urna lub ruletka. Są to w ścisłym znaczeniu tego słowa urządzenia losowe. Generatory takie m ają jednak niewielkie zastosow anie praktyczne i m ogą być przydatne tylko do losow ania niedużych próbek do badań reprezentacyjnych. M ożna zbudow ać tego typu urządzenia w spółpracujące z kom puterem , np. w przeszłości wielokrotnie konstruow ano urządzenia w ykorzystujące zjawisko prom ieniotwórczości lub zjawisko szumów elem entów elektronicznych. Istotnym problem em je st tu jednak problem stabilności takich generatorów: niewielkie zm iany własności fizycznych źródła lub zm iany w arunków otoczenia m ogą pociągnąć za sobą istotne zm iany własności probabilistycznych otrzym ywanych ciągów liczb losowych. W związku z tym generatory fizyczne w ym agają dodatkowych urządzeń testujących i ewentualnie korygujących, co znacznie kom plikuje ich budowę. Pojaw ia się trudny problem synchronizacji okresów testow ania i okresów eksploatacji takich generatorów. W szystkie te kłopoty z jednej strony i łatw ość eksploatacji generatorów program ow ych z drugiej strony spowodowały, że w spółcześnie te ostatnie całkowicie w yparły generatory fizyczne, a pew ne namiastki generatorów fizycznych (np. zegar systemowy) używ ane są tylko do inicjow ania generatora programowego. W szystkie generatory program ow e (ponieważ dalej mówim y tylko o takich generatorach, więc przym iotnik „program owe" będziem y opuszczali), jakie prezentujem y w tym rozdziale, produkują dodatnie liczby całkowite lub bity (liczby O lub 1). N ieujem ne całkowite liczby losowe produkow ane przez rozw ażany generator będziem y oznaczali w zasadzie dużymi literami X ,X 1,X 2,,„, X n ,,„, ale czasami użyjem y innej litery, np. Y lub Z. Te liczby będą zaw sze m niejsze od pewnej ustalonej dodatniej liczby całkowitej m co oczywiście je st związane z arytm etyką 32 komputera: np. w kom puterze 32-bitowym mamy m = 2 . Losowe bity będziem y zwykle

oznaczali przez b 1, b2,.... Liczby z przedziału (0,1), które m ają reprezentow ać liczby losow e o rozkładzie rów nom iernym na przedziale (0,1), będziem y oznaczali literą U (lub V), ewentualnie z odpowiednimi indeksami. O trzym ujem y je zaw sze w w yniku operacji dzielenia U = X/m lub operacji składania bitów losowych w liczbę ułamkową: U =0. bi, b2, . . Za najwcześniejszy algorytm program ow ego generow ania liczb losow ych je st uważany tzw. algorytm kwadratowy von N eum anna (Hammer(1951)). Podstaw ow a idea generatora von N eum anna polega na generow aniu kolejnych N -cyfrow ych (N - parzyste) nieujem nych całkowitych liczb losow ych X n za pom ocą prostej form uły X n = f( X n-I), gdzie funkcja f je st określona w następujący sposób: oblicza się kw adrat liczby X n-I i, ewentualnie dopisując odpow iednią liczbę zer na początku, otrzym uje się w ynik będący liczbą 2N-cyfrową. Z a kolejną liczbę X n przyjm uje się liczbę utw orzoną z N środkowych cyfr tego wyniku. Okazało się jednak, że generator von N eum anna produkuje zbyt krótkie tablice liczb losowych (patrz np. Gajewski i Zieliński (1965)) i z tego powodu został zaniechany. Idea von N eum anna znajduje zastosow anie i rozwinięcie w e w spółcześnie stosowanych generatorach: obecnie używ ane generatory program ow e liczb losowych produkują ciągi liczb X 0,X I,..., przy czym każdy elem ent takiego ciągu je st obliczany za pom ocą ściśle określonej formuły m atematycznej, zastosowanej do pewnej liczby poprzednich elementów. O dnotujm y od razu, że tak tw orzone ciągi liczb m uszą być ciągami okresowymi, co rażąco koliduje z losowością. Oznacza to, że istnieją liczby naturalne v i P takie, że dla i > v mam y X = X i+jp,j = 1,2,.... Fragm ent ciągu X 0,Xi,...;Xv+P- I nazywam y okresem aperiodyczności ciągu, natom iast liczbę P okresem ciągu. W prow adzone pojęcia m ożna zilustrow ać następująco: okres ciągu X 0,X b- •

X v+I, — ,X v+P-I, X v+P^ • •

Okres ciągu zwykle daje się w yznaczyć teoretycznie, chociaż w niektórych Przypadkach m oże to być trudne. M ów im y o tym dokładniej przy szczegółowej prezentacji w ybranych generatorów. Obecnie najczęściej stosowanymi generatorami są: generatory liniowe, generatory oparte na rejestrach przesuwnych, uogólnione generatory Fibonacciego, generatory oparte na odejm owaniu z pożyczką, na m nożeniu 2 przeniesieniem oraz generatory nieliniowe. Przedstaw im y dokładniej w ym ienione klasy generatorów. W tej prezentacji położym y nacisk na podstaw ow e idee konstrukcji, które zilustrujem y przykładami pochodzącym i z najnowszej literatury. Liczba publikacji z tej dziedziny rośnie lawinowo. W śród ostatnio wydanych pozycji przeglądowych poświęconych generatorom liczb losowych o rozkładzie równomiernym chcielibyśm y w yróżnić książkę Tezuki (1995). N ajnow sze inform acje na tem at generatorów (w postaci opisów nowych technik, pakietów program ów w różnych językach program ow ania oraz odnośników do literatury) zaw iera światowa sieć kom puterowa Internet (poprzez serwery W W W , ftp, listy dyskusyjne itp.). Poruszanie się w tym gąszczu inform acji ułatw iają specjalne serwery wyszukujące; jeżeli chodzi o interesującą nas tem atykę generatorów liczb losowych, m ogą to być np. takie słowa kluczowe jak: mndom num ber generators, Simulation, M onte--Carlo.

2.2. Generatory liniowe 2.2.1. Opis Generatory liniow e to generatory postaci Xn+1 = (aXn + ^2Xn-1 + . . . + akXn-k+1 + c) mod m

(2.1)

gdzie do, a 0, a 1, . . . , ak, c i m są ustalonymi liczbami całkowitymi (parametrami generatora), natom iast a mod b oznacza resztę z dzielenia liczby a przez liczbę b. Inicjując działanie generatora, użytkownik dostarcza dane początkowe: X 0,X 1, . . . ,X k. W typow ych im plem entacjach (np. Pascal, C i C++) przyjm uje się k = l, co prowadzi do generatora Xn+1 = (aXn + c) mod m

(2.2)

po raz pierwszy zaproponow anego przez Lehm era (1951). Jeśli c = O, to otrzym ujem y tzw. generator multiplikatywny, a jeśli c ^ O, mówim y o generatorze mieszanym. W dalszym ciągu dokładniej omówim y generatory liniow e postaci (2.2), gdyż z obszernej klasy generatorów liniowych one w łaśnie są standardowo używ ane w e w spółczesnych komputerach. Inform acje na tem at ogólnego generatora liniow ego (2.1) przedstawim y w p. 2.2.4. Zauw ażm y przede wszystkim, że ciągi (2.2) są dość prostymi ciągami determ inistycznym i i wobec tego raczej tylko w w yjątkowych przypadkach m ożem y „udawać", że mam y tu do czynienia z ciągami liczb losowych. Okazuje się, że na drodze czysto arytm etycznych rozważań niektóre z takich ciągów m ożna od razu zdyskwalifikować. Po pierwsze, ciągi produkow ane przez generatory liniow e są ciągami okresowymi. N a przykład, jeżeli w ciągu (2.2) pew na liczba pojawi się po raz drugi (a musi to nastąpić wcześniej lub później, bo istnieje tylko m różnych reszt z dzielenia przez m), to od tej pory cały ciąg będzie już tylko reprodukcją swojego poprzedniego odcinka. Jeżeli okres ciągu je st zbyt mały, to liczby Ci = X /m ,i = 1,2,..., będą zbyt rzadko w ypełniały przedział(0,1) i ciąg takich liczb nie będzie mógł być zaakceptow any jak o ciąg sym ulow anych realizacji zmiennej losowej o rozkładzie równomiernym U(0,1). O tym ja k w ybierać param etry generatora, żeby zagw arantować dostatecznie duży jego okres, pow iem y w p. 2.2.2. Po drugie, okazuje się, że przy ustalonej liczbie d > l punkty (U1, U 2,..., Ud), (U2, U 3,..., Ud+1),... (2.3) ja k również punkty (U1, U 2 , . , Ud), (Ud+1, U d + 2 ,., U2d ) , .

(2.4)

„bardzo nielosowo" w ypełniają kostkę jednostkow ą I d (tzn. przedział [0,1]d w d-wymiarowej przestrzeni R ) . Tu znowu m ożna przez odpowiednie m anipulacje param etram i generatora uzyskać mniej lub bardziej zadow alające wypełnienia; m ów im y o tym w p. 2.2.3. Po trzecie, jeżeli w ynikow y ciąg U1, U2,... m a udaw ać realizację ciągu niezależnych zm iennych losow ych o rozkładzie rów nom iernym U(0,1), to średnia produkow anych przez generator liczb pow inna być rów na 1/2, ich w ariancja 1/12, a w spółczynniki autokorelacji w ciągu tych liczb pow inny być równe zeru. D la ciągów (2.2) znane są teoretyczne w zory opisujące te w ielkości - zależą one od param etrów generatora. Więcej inform acji na ten tem at podam y w p. 2.2.5.

Po czwarte, jeżeli w ynikow y ciąg Uj, U2,... m a udaw ać realizację ciągu niezależnych zm iennych losow ych o rozkładzie równomiernym U(0,1), to ciągi liczb produkowane przez generator powinny spełniać różne testy statystyczne. Sprawie testow ania generatorów poświęcam y cały odrębny rozdział 5. W punkcie 2.2.6 podam y przykładowo kilka generatorów, które pozy­ tywnie przeszły różne kryteria teoretyczne i testy statystyczne.

2.2.2. Okres generatora Okres generatora liniow ego (2.2) je st równy P = min{i: X = X oi > 0} (zauważmy, że param etr aperiodyczności v = 0). Okres ten nie m oże oczywiście przekraczać liczby m. Taki okres przy odpowiednim w yborze param etrów m oże osiągnąć generator mieszany, ale nie m oże go osiągnąć generator m ultiplikatywny. Jeżeli np. m = 2L, to okres generatora m ultiplikatywnego nie przekracza liczby 2l-2, a jeżeli m je st liczbą pierwszą, to m aksym alny okres generatora m ultiplikatywnego je st równy m - 1. Generator m ultiplikatywny ze stałą m = 2L, L > 4, osiąga m aksym alny okres tylko wtedy, gdy X 0 je st liczbą nieparzystą oraz a = 3 m od 8 lub a = 5 mod 8 (zauważmy, że fakt osiągania m aksym alnego okresu nie zależy od liczby L, która tylko decyduje o tym, ja k długi je st ten m aksy­ malny okres). Przykładem generatora o m aksym alnym okresie je st generator RA N D U z param etram i a = 2 1 + 3 i m = 231, który był standardowo używ any w kom puterach IBM 360/370 i PDP11. M a on jednak bardzo krótki okres, a ponadto nie spełnia niektórych testów statystycznych. Innym przykładem je st generator RNB z param etram i a = 22 • 237 + 1 i m = 2 1, opisany, uzasadniony i statystycznie przetestow any w pracy Zielińskiego (1966). Generator m ultiplikatywny z liczbą pierw szą m osiąga m aksym alny okres tylko wtedy, gdy a (m~lj/p ^ 1 mod m dla każdego czynnika pierw szego p liczby m — 1. Przykładem takiego generatora je st generator z param etram i a = 75 oraz m = 231 - l (m je st przykładem tzw. liczby M ersenne'a, czyli liczby postaci 2 - l, g d zie p je st liczbą pierwszą). Generator m ieszany osiąga pełny okres m wtedy, gdy jednocześnie są spełnione trzy następujące warunki: a) liczby c i m nie mają wspólnych dzielników, b) a = l m o d p dla każdego czynnika pierw szego liczby m, c) a = l m od 4, jeżeli 4 je st dzielnikiem liczby m. Przykład takiego generatora uzyskujem y przyjmując: a = 69069, c = l, m = 232. D owody odpow iednich tw ierdzeń m ożna znaleźć w m onografiach Janssona (1966), K nutha (1981) i Ripleya (1987). Zw racam y uw agę jeszcze na jedną w łasność zw iązaną z okresowością ciągów liczb produkow anych przez generatory liniow e w przypadku param etru m nie będącego liczbą pierwszą. M ianowicie, jeżeli liczby X otrzym ywane z takiego generatora zapiszem y w pewnym systemie pozycyjnym, w szczególności w systemie binarnym, w postaci X = b1b2...bL i następnie przez obcięcie początkow ych (najstarszych) bitów utw orzym y nowe liczby, np. X ' = bibi+1...bL, to te nowe liczby także utw orzą ciąg okresowy, a okres nowego ciągu będzie krótszy od okresu ciągu wyjściowego. W szczególności dla c = 0, m = 2L, końcowe (najm łodsze) bity tw orzą ciąg o okresie równym 1; np. jeśli a = 8k+5, X 0 = 4s + l (gdzie k, s są pewnymi liczbami całkowitymi) oraz m = 2 , to trzy końcow e (najm łodsze) bity tw orzą ciąg o wartościach postaci 001 i 101 (Zieliński (1979), A ndersen (1990)). W ynika stąd, że liczby losow e z om awianego generatora liniowego m ogą być używ ane tylko w takich obliczeniach, w których końcow e bity liczb nie odgrywają istotnej roli; w szczególności poszczególne cyfry liczb nie m ogą być traktowane ja k cyfry losowe, chociaż w przypadku „prawdziwych" ciągów niezależnych zm iennych losowych o rozkładzie

rów nom iernym U (0,1) takie postępow anie je st w pełni uzasadnione.

2.2.3. Struktura przestrzenna Jeżeli U i,U 2,... je st ciągiem niezależnych zm iennych losow ych o jednakow ym rozkładzie równom iernym U(0,1), to punkty losow e (2.3) i (2.4) m ają rozkład rów nom ierny na kostkach jednostkow ych ' f w d-wymiarowej przestrzeni R d. Jeżeli natom iast punkty (2.3) i (2.4) są utw orzone przez liczby U 1,U2,... produkow ane przez generator liniowy, to po pierwsze nie w ypełniają one ty c h k o ste k d o sta te c z n ie g ęsto i po d ru g ie - u k ła d a ją się w tych kostkach w regularne struktury geom etryczne. _ Z biór punktów (2.3) je st zaw arty w zbiorze postaci Ld (x + A) gdzie x = (xi, X2 ,...Xd ) xi = 0,

X2 = c / m, Xi = axi-i + c / m , i = 3..., d

oraz A je st kratą, czyli zbiorem w ektorów postaci t 1e 1+ ^ + t ded dla pewnych liniowo niezależnych wektorów bazowych e 1,e 2, . . . , e d i liczb t 1, t 2, . . . , t d przebiegających zbiór liczb całkowitych Z. Przykładow e zbiory Ld (x + A) w raz z w ektoram i bazow ym i, przedstaw iono na rys. 2.1.

Zbiór punktów (2.4) je st oczywiście pewnym podzbiorem zbioru (2.3). Szczególną rolę w opisie rozkładu punktów (2.3) w kostce ' f odgryw a tzw. baza fizyczna {e1,e2,.. . ,ed}, charakteryzująca się tym, że e 1je st najkrótszym wektorem w A oraz dla i > 2 ei, je st najkrótszym wektorem w podprzestrzeni ortogonalnej do podprzestrzeni rozpiętej na wektorach e 1,e2,.. . ,ei-1 (por. rys. 2.1). Jedną z m iar zagęszczenia punktów (2.3) w ' f je st wtedy długość ld najdłuższego w ektora bazy fizycznej, a jed n ą z m iar rów nom ierności rozkładu tych punktów je st iloraz ld/l1, gdzie li je st długością najkrótszego w ektora tej bazy. Inną m iarą gęstości w ypełnienia kostki ' f punktami (2.3) je st m aksym alna odległość D d m iędzy hiperpłaszczyznam i, na których leżą te punkty. Pom ysł obliczania wielkości D d pochodzi z pracy Coveyou i M acPhersona (1967), gdzie zastosow ano m etodę analizy Fouriera, stąd w literaturze problem w yznaczania D d je st nazywany testem spektralnym. Znane są następujące nierówności w iążące w prow adzone w ielkości (Ripley (1987), Tezuka (1995)): l < ld / D d < y d ,

m Dd1yd > 1

da

d >2

gdzie yd je st tzw. s ta łą H e rm ie fa ; a dokładne w artości tej stałej są znane tylko dla d < 8(D2d = 1,4/3,2,4,8,64/3,64,256, d = 1,2,...,8)) (Knuth (1981). O dległość D d m ożna rów now ażnie zdefiniow ać jak o długość najkrótszego niezerow ego w ektora tzw. bazy dualnej do bazy {e1, e2,..., ed}. B azę dualną e * , e*2,...e*d jo k reśla się za

pom ocą zależności: e Te* (gdzie óp je s t d e ltą K ro n e c k e ra : ó^ =1 d la i = j , óp = O d la i f j ) . P o d a n a d e f in ic ja sprow adza obliczanie D d do problem u m inim alizacji całkow itoliczbow ej z k w ad rato w ą fu n k cją celu. W ięcej szczegółów oraz w ybrane algorytm y num eryczne zw iązane z ba daniem struktury geom etrycznej generatorów liniow ych m ożna znaleźć w ob szernej literatu rze (C oveyou i M ac P herson (1967), K nuth (1981), A Serbach i G rothe (1985), F in ck e i P o h st (1985)). Pew ne oszacow ania z góry w ielkości l 1 i Id m ożna uzyskać przez obliczenie długości odpow iednich w ek to ró w dow olnej bazy. T aką bazą m oże być np. b aza z a w ie ra ją c a w ek to r ei = (l, a, a 2,..., a d' 1)/m i w ek to ry e ui > 2, o i-tej składow ej rów nej jed n o ści i pozostałych składow ych rów nych zeru. Lepsze przybliżenie otrzym uje się p opraw iając tę b azę w je d en z n astęp u jących sposobów : (1) W każdej parze w ek to ró w ei,e, w ek to r dłuższy, np. ej, zastępujem y w ektorem e,-sei, gdzie ^ je s t zaokrągleniem liczby eTe /(eTei) do najbliższej liczb y c ałk o w itej, k tó ra je d n o c z e ś n ie b liż s z a je s t zeru. K o n ty n u u je m y tę p ro ced u rę dopóty, dopóki m ożna je szc ze uzyskać redukcję długości rozw ażan y ch w ek to ró w (M arsagli a (1972)). (2) W ektory ei porządkujem y w edług długości i każdy z nich zastępujem y najkrótszym w ektorem postaci e. + / ^ < i c e i gdzie c. e {0,+1,-1}(R ipley (1983)). Stosując pow yższe algorytmy, zaw sze po skończonej liczbie kroków osiągam y bazę, której ju ż nie m ożna popraw ić, gdyż w spółrzędne w ektorów bazow ych są pew nym i w ielokrotnościam i w artości l /m. N a rysunku 2.1 dla dw óch przykładow ych generatorów liniow ych zaznaczono w ektory bazow e zbioru A . W ektory te uzyskano w ykonując je d e n krok zgodnie z algorytm em (1), p rzy czym b aza p o czą tk o w a m iała p o stać: e 1 = (l/m , a/m ), e2 = (0,1).

2.2.4. Ogólne generatory liniowe Jako naturalne u ogólnienie generatorów (2.1) ro zw aża się rów nież ge n e ra to ry p o staci X n+1 = A X n m o d m

(2.5)

gdzie X 1,X 2,... są w ektoram i w R k, k > l, A je s t m acierzą. O peracja m od je st w ykonyw ana „po w spółrzędnych". Te ogólniejsze generatory um ożliw iają sym ulow anie w ielow ym iarow ych zm iennych losow ych o rozkładzie rów no m iernym na kostkach 'f, a m ożliw ość w yboru m acierzy A pozw ala na m anipulow anie zależnościam i m iędzy składow ym i generow anych w ektorów . N ie b ęd ziem y tutaj ro zw ijali tej p ro b lem aty k i; zain tereso w an eg o C zy teln ik a odsyłam y do literatury (L 'E cuye r (1990, 1996a), E ichenauer-H erm ann, G rothe i L ehn (1989)).

2.2.5. Parametry statystyczne Jak już podkreślaliśmy, ciągi liczb z generatora liniowego (2.1) są ciągami deterministycznym i i tylko pew ne ich w łasności „nieuporządkowania" uspraw iedliw iają stosowanie ich do symulacji ciągów losowych: jeżeli w ynikow y ciąg U 1,U2,... m a udaw ać realizację ciągu niezależnych zm iennych losowych o rozkładzie równomiernym U(0,1), to średnia produkow anych przez gene­ rator liczb pow inna być rów na 1/2, ich w ariancja 1/12, a w spółczynniki autokorelacji w ciągu otrzym anych liczb powinny być równe zeru. D la ciągów produkow anych przez generatory liniowe (2.2) znane są w zory teoretyczne dla tych w ielkości (Jansson 1966); zależą one od param etrów generatora. Podam y kilka takich w zorów dla konkretnej klasy generatorów m ultiplikatywnych w celu zilustrow ania zagadnienia. W eźm y pod uw agę dowolny m ultiplikatywny generator liczb L-bitow ych o m aksym alnym okresie,

T2

tzn. o okresie 2 - . Generator taki produkuje tylko cztery rodzaje ciągów: (1) jeśli c = 3 mod 8 o ra z X 0 = 1,3,9 lub 11 mód 16, to cią g X 0, X 1,X 2, ... je st perm utacją liczb postaci 8j + 1 i 8j + 3, j = O, l, 2,..., 2L-3-1; (2) jeśli c = l m od 8 o ra zX 0 = 5,7,13 lub 15mod 16, to c ią g X 0, X 1,X 2,... je st perm utacją liczb postaci 8j + 5 i 8j + 7,_j = O, l,2,..., 2L-3 - 1; (3) jeśli c = 5 mod 8 o ra zX 0=

l m od 4, to c ią g X 0, X 1, X 2,... je st perm utacją liczb postaci

4j + l, j = O, l, 2,..., 2l-2- 1; (4) jeśli c = 5 mod 8 o ra z X 0 = 3 mod 4, to cią g X 0, X 1, X 2, ...je st perm utacją liczb postaci 4j + 3, j = O, l, 2,..., 2l-2- 1. R ozw ażm y np. ciąg postaci (1).Średnia w szystkich liczb U = X / 2 L produkow anych przez generator je st rów na 1/2 - (1/2)L-1, a ich w ariancja 1/12 — 13/(3 • 22m). W ynika stąd, że taki ciąg jako całość wypełnia, przy odpowiednio dużym L, przedział (0,1) prawie tak dobrze ja k ciąg „prawdziwych" liczb losowych. W zory dla w spółczynników autokorelacji są nieco bardziej zawiłe, ale w ykonanie obliczeń w edług tych w zorów nie nastręcza w iększych trudności; np. dla generatora m ultiplikatywnego X n+1 = aX nm o d 2

mamy

y

2“

-1 X V j X V j + ,. -=

- JL (( 2l 2m 2m

6 * 2 m - 1) -

-

,L - 2

2

j =0

^ .i - , j + a /8 ^ „2^ 1/ a j + a /8 - 32 Z j 1 - j L T 4 Z ji 22 L -3 r J j =0 V 2 J j =0 V -

32

-r-n

,|

j=0

i V

Zj

a j

+

3— /8 ^

^L-3r 2

J

. ^ 2^

1i

12

j|

Z

j =0

arj

+

3— /8 L-3 r

V

2

J

r L gdzie ar = a mod 2 . O dpowiednie algorytmy obliczeniowe, a także tablice ju ż obliczonych w spółczynników autokorelacji w standardowych generatorach liniowych, m ożna znaleźć w m onografii Janssona (1966).

Znane są rów nież mniej dokładne, ale w ygodniejsze w projektowaniu generatorów wzory, np. w artość w spółczynnika autokorelacji ^ m iędzy d ^ iem a ^ kolejnymi liczbam i z generatora m ieszanego (2.2) znajduje się w p rz e d z ia le f 1------ | ± — (Greenberger (1961, 1962)).

2.2.6. Wybór parametrów dla generatorów liniowych Ostatecznym kryterium zaakceptow ania generatora je st to, że nie został on zakw estionowany przez żaden z zastosow anych testów statystycznych. M etody testow ania om ówim y dokładnie w rozdz. 5., gdzie również opiszem y liczne testy statystyczne. Tutaj, mówiąc o generatorach liniowych, zwracam y uwagę, że jeżeli wartości pew nych parametrów, ja k np. średnia i w ariancja produkow anych liczb, okres ciągu lub w spółczynniki autokorelacji w ciągu

liczb produkow anych przez generator, m ogą dla danego generatora być w yznaczone teoretycznie w sposób ścisły, to oczywiście nie m a żadnego sensu testowanie hipotez o takich param etrach za pom ocą testów statystycznych. W obszernej, liczącej ju ż ponad 40 lat, literaturze poświęconej generatorom m ultiplikatywnym i mieszanym, m ożna znaleźć liczne raporty z wykonanych analiz teoretycznych oraz testów statystycznych i w yróżnić generatory, które najlepiej spełniały przyjęte kryteria porównawcze. W tabeli 2.1 przedstawiam y kilka takich generatorów postaci (2.2); odnośniki do w ielu w cześniejszych prac znaleźć m ożna w książce Zielińskiego (1979). W szytkie podane w tabeli generatory osiągają m aksym alne okresy. Popularność param etrów m postaci 2 2 lub 231 — l w ynika z łatw ości im plem entacji takich generatorów w e w spółczesnych systemach komputerowych. W iele generatorów z tabeli nadal stanowi standardowe wyposażenie szeroko używ anych system ów operacyjnych, języków programowania, a naw et specjalistycznych pakietów oprogram ow ania do obliczeń naukowych. N ależy tutaj podkreślić, że w obecnie realizow anych obliczeniach sym ulacyjnych z użyciem liczb pseudolosow ych okres generatora rzędu 232 je st zbyt mały, gdyż zużycie w szystkich liczb z takiego generatora następuje za szybko. W ielu autorów sugeruje, aby liczba N liczb z generatora używ anych w symulacji była dużo m niejsza niż okres generatora P. W pracy M aclarena (1992) uzasadniano, że liczba N nie pow inna przekraczać liczby P 2/3; odpowiednie ograniczenie w przypadku generatorów liniowych, podane przez Ripleya (1987), wynosi P 1/2 Generatory liniow e postaci (2.2) nie spełniają pew nych now szych testów statystycznych, np. testu OPSO (patrz rozdz. 5.). Eksperym enty num eryczne potwierdziły rów nież (np. M arsaglia (1984,1995)) lepsze w łasności statystyczne generatorów liniow ych konstruow anych z param etrem m będącym liczbą pierw szą (jednak kom plikuje to im plem entację generatora oraz w pływ a na jego szybkość).

T ab e la 2.1

a

c

m

Źródło

22 • 237 + 1

0

235

Zieliński (1966)

69069

1

232

M arsaglia (1972)

16807

0

231-1

Park, M iller (1980) Carta (1990)

630360016 397204094

0

2 31-1

0

231-1

410092949

0

232

Borosh, N iederreiter (1983)

742938285

0

231-1

Fishman, M oore (1986)

40692

0

231 - 249

L'Ecuyer (1988)

1099087573

0

232

Fishm an (1990)

68909602460261

0

248

Fishm an (1990)

Fishman, M oore (1982)

Z ogólniejszych generatorów liniow ych (2.1) dobrą ocenę statystyczną (w tym w e wspom nianych now szych testach) uzyskały m.in. następujące generatory (M arsaglia (1995)):

1) X n = (1176Xn-1 + 1476Xn_i+ l776Xn-j; m od (232 - 5) 2) Xn = 2 13(Xn-1 + Xn-2 + Xn-3) mod (232 - 5) 3) X n = (1995Xn-1 + 1998Xn-2+ 2001Xn-3) m od (235 - 849) 4) X n = 2 19(Xn-1 + Xn-2 + Xn-3) mod (235 - 1629) Generatory te osiągają m aksym alne okresy równe m — l, gdzie liczba m je st odpowiednim 32 modułem: dla pierw szego i drugiego z pow yższych generatorów równym 2 -5, dla trzeciego 235-849 oraz dla czwartego 23 -1629. Przedstaw ione generatory m ogą być łatw o zaim plem entowane na w spółczesnych kom puterach w każdym języku programowania. Przykładow ą im plem entację ogólnego generatora liniow ego zam ieszczam y w podrozdz. 2.11 (patrz rów nież uwagi w podrozdz. 2.10).

2.3. Generatory oparte na rejestrach przesuwnych P rzyjm ijm y, że k je s t ustalo n ą liczbą n aturalną i w eźm y pod uw agę ciąg bitów zdefiniow any w zorem rekurencyjnym bi = (aibi- 1 + ... + a kbi-k)m od2,

i = k + 1,k + 2,...

(2.6)

gdzie w spółczy n n ik i a 1,a2..., ak są stałym i b in arn y m i, tzn. liczbam i 0 lub l, oraz b 1,b2, . .., b k je s t ustalonym ciągiem inicjującym . Z ależność (2.6) m ożna opisać za pom ocą operatora logicznego xor, zw anego różnicą sym etryczną lub alternatyw ą w yłączającą, o następującej tabelce działań: a

b

a xor b

0 0 1 1

0 1 0 1

0 1 1 0

Inaczej to ujmując, dla zmiennych boolowskich mamy po prostu a xor 6 = (a + 6) mód 2 Dalej będziemy używać tego operatora również w odniesieniu do liczb całko witych; w takim przypadku wynik jest liczbą całkowitą złożoną z bitów powstałych w wyniku działania operatora na poszczególnych pozycjach reprezentacji binarnej argumentów. Wzór (2.6) przyjmuje wtedy równoważną postać: jeśli aji =...=aj, = l, a pozostałe współczynniki są równe zeru, to bi = bi-ji xor b j xor ... xor bH, Ciąg bitów (2.6) jest oczywiście ciągiem okresowym. Ponieważ istnieje 2k różnych układów kelementowych (b1,b 2, ... ,bk), okres ciągu (2.6) nie może być większy od 2k. Faktycznie może on być równy co najwyżej 2k - l, bo gdyby pojawiło się w nim k kolejnych zer, to cały ciąg musiałby składać się z samych zer. Mówiąc dalej o ciągu o maksymalnym okresie, mamy na myśli ciąg o okresie 2k - 1. Metodami algebraicznymi można badać, dla jakich współczynników ai, a 2, ..., ak ciąg (2.6) ma maksymalny okres. Nie będziemy tutaj rozwijali tego technicznego wątku (odsyłamy Czytelnika np. do pracy Golomba (1967)), ograniczymy się natomiast do prostszego, ale dla nas atrakcyjnego przy padku, gdy w formule rekurencyjnej (2.6) tylko dwa współczynniki at są różne od zera. Schemat iteracyjny (2.6) ma wtedy postać bi = b-p xor bi-q

(2.7)

dla pewnych ustalonych liczb naturalnych p oraz q. Bez zm niejszania ogólności rozważań przyjmijmy, ż e p > q. W tabeli 2.2 podajem y przykładowe wartości param etrów p i q, dla których ciąg (2.7) ma m aksym alny okres. Tabelę tę utw orzono na podstawie tabeli z książki Ripleya (1987), gdzie podano wartości (p,q) zapewniające m aksym alny okres d la p < 36 oraz tabeli z Berdnikova, Turtii i Com pagnera (1996), w której z kolei zawarto znane wartości (p,q), dające maksym alny T ab e la 2.2

P

q

P

q

2

1

33

13

3

1

35

2

4

1

36

11

5

2

89

38

6

1

127

1, 7, 15, 30, 63

7

1,3

521

32, 48, 158, 168

9

4

607

105, 147, 273

10

3

1279

216, 418

11

2

2281

715, 915, 1029

15

1,4,7

3217

67, 576

17

3,5,6

4423

271, 369, 370, 649,

18

7

20

3

9689

84, 471, 1836, 2444, 4187

21

2

19937

881, 7083, 9842

22

1

23209

1530, 6619, 9739

23

5,9

44497

8575, 1034

25

3,7

110503

25230, 53719

28

3, 9, 13

132049

7000, 33912, 41469,

29

2

31

3, 6, 7, 13

1393, 1419, 2098

52549, 54454

okres ciągu (2.7) dla liczb p < 132049, gdzie 2P - l są liczbami pierwszymi (założenie to znacznie upraszcza obliczenia). D odatkow o m ożna rozważać pary w spółczynników (p,p — q), gdyż jeśli ciąg (2.7) m a maksym alny okres dla pew nychp i q, to własność tę m a rów nież schemat iteracyjny z param etram i p i p — q. Istnieje wiele sposobów uzyskiw ania L-bitow ych liczb losowych o w artościach w przedziale (0,1) na podstaw ie ciągu bitów (bj,i = 1,2,...). N ajprostszy polega na konstruow aniu ich za pom ocą wzoru 2. G eneratory o rozkładzie równomiernym L U i = ^ 2 - bis+j = 0 bis+1 ■■■bis+L, j =1

i=0,1,2,...

(2.8)

gdzie ^ je st ustaloną dodatnią liczbą całkowitą oraz s < L. Jeżeli s < L, to kolejne liczby Ui w ykorzystują te same podciągi bitów, a jeżeli s — L, to liczby Ui i U i+1 utw orzone są z rozłącznych fragm entów ciągu (bi = 1,2,...). Liczby U i otrzym uje się łatw o za pom ocą rejestrów przesuwnych oraz bram ek logicznych, realizujących operator xor (ten fakt uzasadnia nazw ę rozważanej klasy generatorów). Generator (2.8) je st znany w literaturze jak o generator Tauswortha, gdyż po raz pierwszy był analizow any w pracy Tausw ortha (1965). Jeżeli liczba s je st tak wybrana, że nie ma w spólnych dzielników z liczbą 2k - l, to ciąg (2.8) je st ciągiem o m aksym alnym okresie, który jest równy 2k - 1. Efektywny algorytm generow ania tak zdefiniow anych ciągów (Ui, i= 1,2,...) za pom ocą ciągu bitów (2.7), w przypadku, gdy q < p/2 oraz O < s < p - q, podał w swojej m onografii Tezuka (1995). W formalnym zapisie tego algorytmu użyjem y następujących symboli: A > k oznacza takie przesunięcie w prawo. Przy przesuw aniu w lewo m łodsze (zwalniane) bity są zapełniane zerami. Przy przesuw aniu w prawo zw alniane starsze bity są rów nież zapełniane zerami. Zw racam y jednak uw agę na to, że w konkretnej im plem entacji działanie operatora przesunięcia bitowego w prawo m oże zależeć od zadeklarow anego typu liczby całkowitej.

W poniższym algorytm ie A je st L-bitow ą liczbą całkowitą o bitach początkow ych bi,..., &l (są to bity składające się na liczbę C/o), natom iast B je st L-bitow ą zm ienną pomocniczą. A LG ORY TM T:IM PLEM EN TA CJA SCHEM ATU TA U SW ORTHA

(2.8)

D LA 0 < s < p -q 1: B = ((A « q) xor A) « (L - p) 2: A= (A« s) xor (B » (L - s)) 3: Return A; goto l

Algorytm ten je st bardzo łatw y w realizacji np. w języku C, który zaw iera operatory przesuw ania bitowego oraz operator xor. W podrozdziale 2.11 zam ieszczam y przykładową im plem entację w tym języku. Podana tam realizacja używ a dodatkowo kom binacji trzech generatorów om awianego typu za pom ocą operatora xor (patrz podrozdz. 2.6). W książce Tezuki (1995) je st zaw arta analiza struktury punktów w kostce (O,1)d, d = 2,3,...,20, tw orzonych z kolejnych liczb z generatora. Lew is i Payne (1973) zaproponowali inny schemat generow ania L-bitow ych liczb całkowitych Yi na podstawie ciągu (2.6), w edług zależności Yi = bibi-l2...bi-Ł

gdzie l2, . . .lL są ustalonymi param etram i przesunięcia. D la ciągu (2.7) otrzym ujem y stąd schemat Yi = Y i-p xor Yi-q. Generatory realizujące ten schemat są nazyw ane uogólnionym i generatoram i opartym i na rejestrach przesuwnych. W ym agają one odpowiedniego w yboru punktów startowych (Y0,Yi, ..,YP) (Lewis i Payne (1973), Collings i H em bree (1986)). Interesujące przykłady efektywnej im plem entacji w omawianej klasie generatorów zaw ierają artykuły: K irkpatricka i Stolla (1981) (dla p = 250, q = 147) oraz Ripleya (1990) (dla p = 521, q = 32). W pracy Berdnikova, Com pagnera i Turtii (1996) zaproponow ano kom binacje (patrz podrozdz. 2.5) kilku generatorów z omawianej klasy, dla odpowiednio dobranych param etrów p i q. G warantuje to bardzo duży okres generatora (np. dla konkretnej kombinacji czterech generatorów 16378 składowych uzyskano okres rzędu 10 ). Ponadto wykazano, że taka konstrukcja zapew nia

bardzo dobre własności dotyczące niezależności kolejnych długich (rzędu kilkunastu tysięcy) ciągów liczb z generatora (patrz rów nież Com -pagner i W ang (1993), Com pagner (1995)). Kody źródłow e w języku C om awianych generatorów m ożna znaleźć w Internecie (np. http://www.can.nl lub ftp.can.nl) .

2.4. Generatory Fibonacciego Ciąg rekurencyjny f n = f n-2 + f n-\-, n > 2 f 0 = f 1= 1 badał ju ż Fibonacci (Leonardo z Pizy) i wyniki swoich badań opublikował w pracy Liber abaci w 1202 roku. Ciąg reszt, przy ustalonej dodatniej liczbie całkowitej m, zdefiniow any wzorem X n = X n.2 + X nl m o d m,

n > 2,

(2.9)

zachowuje się na tyle bezładnie, że ju ż dawno zainteresował m atem atyków poszukujących prostych modeli dla procesów losowych. W ydaje się, że pierw sze wyniki sprawdzania tego ciągu za pom ocą testów statystycznych zostały opublikowane w pracy Taussky i Todd (1956). Okazało się, że ciągi (2.9) spełniają testy równom ierności rozkładu, ale nie spełniają testów niezależności, a w ięc także w ielu innych testów (np. testów serii), gdzie niezależność odgrywa kluczow ą rolę. Tej w ady m ożna było się pozbyć, uogólniając ciąg (2.9): X n =• X n_r + X n-8 m od m ,

n>r,

r>s> 1

(2.10)

ale odbywało się to kosztem czasu, co czyniło takie generatory mało konkurencyjnymi dla rozpow szechnionych generatorów multiplikatywnych. Argument, że takie uogólnione generatory miały dużo dłuższy okres od generatorów m ultiplikatywnych, nie był przekonujący, bo na stosunkowo w olnych kom puterach rzadko dochodziło do w yczerpania okresu ciągów m ultiplikaty­ wnych. Trudno się natom iast dziwić, że w dzisiejszej dobie szybkich kom puterów generatory Fibonacciego, w różnych w ersjach uogólnień, przeżywają praw dziw y renesans. N astępny krok w tych uogólnieniach polega na zastąpieniu dodaw ania w (2.10) jakąś inną operacją. Ogólnie, w ybraną operację oznaczam y symbolem o i zakładamy, że je st ona wykonywana m odulo m. Uogólniony generator Fibonacciego przyjm uje w tedy postać X n = X n-r ◊ X n-s, i je st oznaczany przez F(r, s,0).

n > r,

r > s >l

(2.11)

r LI , to m aksym alny okres generatorów F(r,s, +) oraz F(r,s,-) je st równy (2 - l)2 - ; dla F (r,s,*), gdzie * oznacza tu (i w całej książce) mnożenie, wynosi on (2r - 1)2n-3, natom iast dla generatora F(r,s, xor) wynosi 2r -1. D owodzi się tego, korzystając z pewnych własności odpow ied­ nich macierzy; szczegóły m ożna znaleźć w pracach M arsaglii (1984) oraz M arsaglii i Tsaya (1985). Generatory F(r,s,xor) dla m = 2 opisaliśmy ju ż w poprzednim podrozdziale, w klasie generatorów opartych na rejestrach przesuwnych. Poniżej podajemy tabelkę z przykładowym i param etram i, zapewniającym i maksym alny okres generatora (2.11): /

r 17 31 55 68 97 607 1279

5 13 24 33 33 273 418

N a przykład, generatory F(17,5,+) oraz F(17,5,-) dla m = 232 dają ciągi o okresie (217 - 1)231, generator F(17,5,*) osiąga okres (2 17 - 1)229, natom iast F(17,5, xor) osiąga okres 2 17-1 = 131071. Omawiane generatory są łatw e w implementacji.

2.5. Kombinacje generatorów D ośw iadczenia z użyciem generatorów skonstruowanych przez łączenie dwóch lub większej liczby prostszych generatorów wykazały, że generatory takie m ają lepsze własności statystyczne niż generatory wyjściowe. Przytoczym y wyniki, które potw ierdzają te obserwacje. Załóżmy, że mamy zm ienne lo so w eX i Y, określone na zbiorze S = {1,2,..., n}, o rozkładach praw dopodobieństw a P (X = ij= p l,

P { Y = i} = ql, i =1,2,...,n

N iech p = l,... ,ro będzie dowolną ustaloną liczbą i niech t = (t1,t2,...,tn) będzie ustalonym wektorem. W eźm y pod uw agę p-norm ę tego wektora, zdefiniow aną wzorem ( n Y 'P IU II= |l t1p v i=1 y

D la zdefiniowanej wyżej zmiennej losowej X w prow adźm y m iarę ó(X) „bliskości" rozkładu tej zmiennej do rozkładu rów nom iernego na zbiorze S Ó(X) = ||(pl,p2,... ,pn) - (l/n, l/n,..., l/n)|| Rozpatrzm y dwuargum entowe działanie o na zbiorze S takie, którego tabelka tw orzy kw adrat łaciński (tzn. każdy jej w iersz i kolum na je st pew ną perm utacją elem entów zbioru S). M ożna udow odnić (np. Brow n i Solomon (1979)), że rozkład zmiennej losowej X o Y je st bliższy rozkładowi rów nom iernem u na zbiorze S w tym sensie, że Ó(X) < m in {Ó(X),Ó(Y)} W odniesieniu do ciągów produkow anych przez generatory liczb losow ych ten w ynik teoretyczny m ożna interpretow ać w następujący sposób: jeśli tablica działań operatora o je st kwadratem łacińskim , to nowy ciąg (X1 o Y1, X 2 o Y2,...) pow inien być bardziej rów nom iernie (a przynajmniej nie mniej równom iernie) rozłożony niż każdy z ciągów X 1, X 2,...i Y1, Y2,... Najczęściej za operator o przyjm uje się rozw ażane wcześniej operatory +, -, *, xor. Ponadto okazuje się, że kombinacje

generatorów produkują ciągi, które są nie tylko „bardziej równomierne", ale rów nież „bardziej niezależne". Podaną konstrukcję dla dwóch generatorów w naturalny sposób uogólnia się na w iększą liczbę generatorów składowych. Ponadto kom binacje generatorów produkują ciągi o w iększym okresie niż okresy ciągów składowych. W szczególności wiadom o, że jeśli ciąg X 1, X 2,... m a okres P 1 oraz ciąg Y1,Y2, . m a okres P 2, gdzie P 1 i P 2 są liczbam i względnie pierwszymi, to okres ciągu X 1oY1,X2 o Y2,... wynosi P 1P 2. Fakt ten je st prostym w nioskiem z tzw. chińskiego twierdzenia o resztach (patrz np. Graham, Knuth, Patashnik (1996)).

Idea kom binow ania generatorów w celu zw iększenia okresu i polepszenia własności statystycznych m a ju ż swoją historię za sobą, ale ciągle je st bardzo często stosowana w najnowszych konstrukcjach generatorów. Pierw sze pom ysły dotyczące kombinacji generatorów pojawiły się w pracach M ac Larena i M arsaglii (1965) oraz M arsaglii i B raya (1968). D użą popularnością w śród użytkow ników cieszył się pomysł kom binowania generatorów za pom ocą dodaw ania m odulo l w dziedzinie liczb rzeczywistych, zaproponow any w pracy W ichm anna i H illa (1982). W yniki teoretyczne związane z kombinacjami generatorów m ożna znaleźć rów nież w pracach: Brow n i Solomon (1979), M arsaglia (1984), Deng (1990), L'Ecuyer i Tezuka (1991). N ajnow sze rezultaty zaw iera np. praca L'Ecuyera (1996b).

2.6. Uniwersalny generator liczb losowych o rozkładzie równomiernym Przez generator uniwersalny rozum iem y generator dający identyczne wyniki na komputerach, w których liczby całkowite są reprezentowane, przez co najmniej 16 bitów, a liczby w arytm etyce zm iennopozycyjnej m ają przynajmniej 24-bitową reprezentację mantysy. Przedstaw im y tutaj szczegółowo generator liczb losowych o rozkładzie rów nom iernym na przedziale [0,1), opublikowany w pracy M arsaglii, Zam ana i Tsanga (1990). Będziem y go nazywali generatorem M ZT. Spełnia on w szystkie znane testy statystyczne i m a duży okres, rów ny 2 144. Generator ten je st kom binacją dwóch prostszych generatorów. Pierw szy z nich je st generatorem typu F(97, 33, o) (patrz podrozdz. 2.4) i produkuje liczby Vn z przedziału [0,1) w edług w zoru rekurencyjnego: Vn = Vn-97 o Vn-33 gdzie x o y = x - y dla x > y oraz x o y = x-y + l dla x < y. Zainicjow anie generatora polega na w yznaczeniu liczb V 1, V2,... ,V97 za pom ocą ciągu bitów (bn) w taki sposób, że V 1= 0.b1b2... b24, V2 = 0.b25b26...b 48, ... itd. Z kolei ciąg bitów je st generow any za pom ocą kombinacji dwóch różnych i łatw ych w im plem entacji generatorów zadanych ciągami liczb całkowitych (yn) i (zr): yn= (yn-3*yn-2*yn-l) m od 179 z n = (52zn1+1) mod 169

0, bn = 1

je ś liy nz n mod 64 < 32 w przeciwnym w ypadku

W ybór m ałych wartości 179 i 169 zapewnia uniw ersalność generatora. U żytkow nik musi dostarczyć procedurze inicjow ania czterech w artości całkowitych:

y 1y 2y 3 € {l, 2,..., 178} (nie w szystkie równe 1) , z 1 G {O, l,..., 168} 120 Okres tego generatora je st równy 2 . Drugim generatorem je st generator liczb losow ych z przedziału (0,1): cn = cn1 o (7654321/16777216),

n > 2,

cl = 362436/16777216

gdzie c o d = c - d dla c > d oraz c o d = c - d + (16777213/16777216) dla c < d, c, d G [0, 1). Okres tego generatora je st rów ny 224 - 3. O statecznie generator M ZT przyjm uje postać Un

Vn o cn

Przykładow ą im plem entację (w języku C) generatora M ZT podam y w podrozdz. 2.11.

2.7. Generatory oparte na odejmowaniu z pożyczką i generator ULTRA W pracy M arsaglii i Zam ana (1991) w prow adzono now ą klasę generatorów wykorzystującą operację tzw. odejmowania z pożyczką (SWB - ang. substract w ith borrow). W tej operacji (oznaczym y ją przez 0 ) bierze dodatkowo udział param etr c, przyjm ujący wartości O lub l, zwany bitem przeniesienia. W ynikiem operacji x 0 y m od m są liczby x -y - c + m lub

oraz

x - y - c oraz

c=l, c=O

gdy x - y - c < O w przeciwnym przypadku,

a początkow ą w artością bitu przeniesienia c je st 0.

Omaw iana klasa generatorów została w ykorzystana przez M arsaglię i Zam ana do konstrukcji generatora ULTRA. Parametram i generatora U LTRA są dodatnie liczby całkowite L > 32, s, r (r > s). W ygodnie je st używ ać oznaczenia m = 2 . W etapie inicjow ania generatora tw orzy się ciąg liczb całkowitych X 1, .

X

e (O, m) oraz ustala się c = 0. W etapie roboczym

w yznacza się kolejne liczby X n w edług w zoru

X n = X n-r 0 X n-s mod m Etap inicjow ania polega na utw orzeniu L-bitow ych liczb X 1, X 2,... ,Xr: X i = bibL-1.b1

(2.12)

X2 = b2lb2L-1...bL+1 Gdzie b1,b2, .

je st ciągiem utw orzonym z bitów znaków liczb w 1,w 2,..., otrzym ywanych w

następujący sposób: użytkow nik generatora podaje dwie liczby u0,vo e (O,m), a następnie są obliczane kolejno liczby:

u1 = rk ui-1 mod m Vi = V >> k 1)x o r Vi Vi = V > 16); y = 30903 * (y&65535) + (y >> 16); return ((x < < 16) + (y&65535)); gdzie zm ienne x oraz y oznaczają liczby 32-bitowe, zawierające w swych starszych i m łodszych 16bitowych częściach odpow iednio wartości X n,c1 oraz Yn,c2, natom iast & je st standardowym operatorem języka C.

2.9. Generatory nieliniowe Przedstaw ione w poprzednich podrozdziałach klasy generatorów są oparte na liniowych wzorach rekurencyjnych. N iepożądaną konsekw encją tej liniowości je st fakt, że odpowiednie punkty (2.3) i (2.4) w przestrzeni wielowym iarowej skupiają się tylko na pewnej liczbie hiperpłaszczyzn, co rażąco odbiega od naszych oczekiwań wobec punktów losowych (por. p. 2.2.3).

N a przezwyciężenie przeszkód naturalnym pom ysłem w ydaje się zatem rozw ażenie ciągów opartych na formułach, które nie są liniowe. Ten kierunek badań nad generatorami rozw ija się dopiero od niedaw na i nadal je st otwarty. Przedstaw im y tutaj pew ne wyniki, które dotyczą generatorów opartych na obliczaniu odwrotności oraz kwadratów. W pracy Eichenauera i Lehna (1986) zaproponow ano generator X n+1 = (aXn 1 + b) mod m,

n = 0,1,...

(2.13)

gdzie m je st liczbą pierwszą. Odwrotność m odulo m je st definiowana następująco: jeśli c = 0, to c-1 mod m = 0, w przeciwnym przypadku c '1 mod m je st taką liczbą całkowitą, że c*c'1 mod m = l, czyli c '1 = cm-2mod m. W ten sposób uzyskujem y ciąg o w artościach w zbiorze {O, l,... ,m - 1}, który m ożna w zwykły sposób przekształcić na ciąg liczb w przedziale [0,1) za pom ocą w zoru Un = X n/m. Innym w ariantem generatora opartego na odwrotności je st zaproponow any w pracy Eichenauera-H erm anna (1993a) generator X n = (a(n + n0) +b) -1 m od m,

n = 0,1,...

(2.14)

który w yróżnia się tym, że kolejna w artość X n może być uzyskana niezależnie od innych elem entów produkow anego ciągu. Taki generator może być szczególnie użyteczny w obliczeniach prow adzonych na kom puterach równoległych. O kazuje się, że dla każdej liczby a e {1,2,..., m) jego okres je st równy m, a w ięc je st to generator o okresie maksymalnym. W arunek na to, aby generator (2.13) osiągał m aksym alny okres równy m, je st nieco bardziej skomplikowany: tak się dzieje wtedy, gdy m 2 - l je st najm niejszą liczbą całkowitą taką, że zm2"1 - l mod (z2 - bz -a). Struktury przestrzenne tw orzone przez odpow iednie ciągi w przestrzeniach wielowym iarow ych dla omawianych generatorów są opisane przez Eichenauera-H erm anna (1991) oraz N iederreitera (1994)). Ponadto w iele eksperym entów obliczeniowych potw ierdziło dobre własności statystyczne tej nowej klasy generatorów. W pracach Eichenauera- H erm anna (1993b, 1994, 1995) analizuje się kom binacje kilku generatorów postaci (2.13) lub (2.14) (patrz rów nież podrozdz. 2.5). Przypuśćmy, że mam y r > 5 takich generatorów Xn(i),

j = 1,... ,r

z param etram i m , j = l , . , r, które są liczbam i pierwszymi. Zdefiniujm y nowy ciąg Un= (Un0 + ...+ Un(r)) m od 1,

n = 0,1,2,...

gdzie Un0 = X n(l^>/mi, j= l, 2 , . , r. Okres tego ciągu je st rów ny m = m f ...^ m r. Pozw ala to na efektywną konstrukcję generatora o długim okresie, przy czym - ja k się okazuje - ten wynikowy generator zachowuje dobre własności strukturalne generatorów składowych. Innym przykładem generatora nieliniowego je st rozw ażany w pracy Blum a (L.), Blum a (M.) i Shuba (1986) generator liczb losowych postaci X n+1 = X n mod m ,

n = O, l, 2 ...

W cytowanej pracy pokazano zastosow ania takiego generatora w kryptologii (jest to dziedzina badań zajm ująca się m etodami szyfrowania informacji). Bardziej szczegółowe inform acje

dotyczące zw iązków m iędzy teorią szyfrowania a generatorami liczb losowych, Czytelnik może znaleźć na przykład w książkach Tezuki (1995) oraz Schneiera (1995).

2.10. Uwagi o implementacji numerycznej Jak już m ów iliśm y (p. 2.2.1), generatory produkują ostatecznie liczby U1,U 2, . z przedziału (0, 1). W niektórych generatorach m ogą pojawić się kłopoty, gdy w ciągu kolejno otrzym ywanych liczb pojawi się O lub 1. Te kłopoty m ogą mieć charakter w ew nętrzny w tym sensie, że po w yprodukow aniu zera generator w dalszym ciągu produkuje ju ż tylko same zera, albo zewnętrzny w tym sensie, że niektóre działania na kolejnych liczbach losow ych m ogą okazać się niewykonalne (dzielenie przez zero lub liczbę bliską zeru, logarytm ow anie takiej liczby, itp.). D odatkow a trudność polega na tym, że faktycznie liczby U1,U2,... są uzyskiw ane z liczb całkowitych X1,X2,... za pom ocą standardowej operacji dzielenia Un = X n/m i w ówczas m ała liczba całkowita X n m oże zam ienić się w m aszynowe zero. Zw racam y uw agę na ten rodzaj trudności, chociaż nie potrafim y podać tutaj żadnej uniwersalnej recepty na poradzenie sobie z nimi. Inny problem je st zw iązany z tym, że typow e generatory są oparte na arytm etyce reszt względem ustalonej dodatniej liczby całkowitej m. Przed obliczeniem takiej reszty m usim y jednak wykonać pewne inne operacje, np. dodawanie lub mnożenie, na liczbach całkowitych ze zbioru {1,2, ...,m), a otrzym any w ynik może w yprowadzać i z reguły w yprow adza poza ten zbiór. Istnieją specjalne algorytm y pozw alające operować liczbam i całkowitymi w taki sposób, żeby wyniki działań pośrednich nie przekraczały liczby m. N ie będziem y rozw ijać dalej tego wątku; zainteresow anego Czytelnika odsyłam y do prac: L'Ecuyer (1990), L'Ecuyer i Cote (1991), D wyer (1995), K nuth (1981).

2.11. Przykładowe implementacje w języku C 2.11.1.

Inicjowanie generatorów

W konkretnych realizacjach algorytm ów generujących liczby losowe, pojawia się problem w yboru odpowiednich wartości inicjujących generator. Jeżeli użytkownik chce zautom atyzow ać rów nież ten początkow y etap obliczeń, może wykorzystać jakiś generator fizyczny np. zegar systemowy. W standardowym języku C wbudowany generator liniow y inicjuje się np. w ywołując funkcję srand((u n sig n ed )tim e(N U L L )); liczba początkow a dla generatora je st wtedy równa liczbie sekund, które upłynęły od l stycznia 1970 roku. „Bardziej losowe" liczby startowe dla generatora m ożna uzyskiw ać za pom ocą dostępnych w systemie param etrów związanych z czasem oraz datą. A nderson (1990) zaproponow ał następującą form ułę dla liczby początkowej X 0: X 0 = r + 100(m - l + 12(d - l + 31(g + 24 * (min + 60 * s)))) gdzie r oznacza dwie ostatnie cyfry roku, m - m iesiąc (od l do 12), d - dzień (od l do 31), g godzinę (od O do 23), m in - m inutę (od O do 59) oraz s - sekundę (od O do 59). Pew ną odm ianą tej formuły je st podana w opracow aniu Rundom num bers (1991-1995) propozycja X 0 = s + 60(min + 60(g + 24(d - l + 31 (m - l + 12r)))) Dodatkowo, w celu zapew nienia nieparzystości liczby X 0, zaleca się zam ianę ostatniego bitu tej

liczby na 1. Powyższe m etody są z reguły stosowane do prostych generatorów liniowych, natom iast w przypadku innych klas generatorów, wym agających dużej liczby punktów startowych, zaleca się inicjow anie za pom ocą prostszego generatora, np. liniowego, który zainicjowano opisanymi wyżej metodami.

2.11.2.

Ogólny generator liniowy

Om ówiony tu generator je st im plem entacją schematu liniow ego (2.1). Okres takiego generatora jest rzędu 296. Ponadto generator ten spełnia znane testy statystyczne, m.in. baterię testów D IEHARD (M arsaglia (1995)). W proponowanej im plem entacji plik nosi nazw ę ecng.c. D o inicjow ania ge­ neratora służy procedura init_ecng(), która w ykorzystuje trzy nieujem ne liczby całkowite i, j, k podane przez użytkownika. G enerator je st realizow any przez funkcję ecng(), która przed pierwszym użyciem w ym aga inicjow ania za pom ocą procedury init.ecng(), a której w ynikiem je st liczba losow a typu double (jest to przekształcona na typ double i przedział (0,1) liczba całkowita 32-bitowa). O dpowiednia im plem entacja przyjm uje następującą postać:

/* PLIK: ecng.c */ »include /* zmienne wykorzystywane przez generator */ static double a,b,c; void init_ecng(int ia, int ib, int ic) { a = ia; b=ib; c=ic; } double ecngO { static int n; static double d,x; d = (a + b + c) * 8192; x =fmod(d,4294967291.0); a = b; b=c; c=x; if (x < (float)2147483648.) n=(int) x; else n = (int) (x - 4294967296.); return (n*2.3283064365e-10) ; {

2.11.3. Generator Tauswortha z podrozdziału 2.3 Omawiany tu generator je st im plem entacją algorytmu opisanego w podrozdz. 2.3, a dokładniej je st to kom binacja za pom ocą operatora xor trzech generatorów opartych na rejestrach 28 29 31 26 przesuwnych. Okres generatora wynosi (2 -1) (2 -1) (2 -1), czyli je st rzędu 3 • 10 . Pewne własności teoretyczne takiej konstrukcji zostały podane w książce Tezuki (1995). Proponowany generator spełnia znane testy statystyczne. W prezentowanej im plem entacji plik nosi nazwę tezuka.c. D o inicjow ania generatora służy procedura init(), która w ykorzystuje trzy nielokalne nieujem ne liczby całkowite podaw ane przez użytkownika: sl, s2, s3. Liczby te pow inny spełniać 28 29 31 nierówności: s l < 2 , s2 < 2 , s3 < 2 . Generator realizow any przez funkcję com bT(), która przed pierwszym użyciem w ym aga zainicjowania za pom ocą procedury init(), a której wynikiem jest liczba losow a typu double (dokładniej m ówiąc precyzja w yniku wynosi 32 bity), przyjm uje postać:

/* PLIK: tezuka.c */ /* zmienne wykorzystywane przez generator */ static unsigned int s1, s2, s3; void init(unsigned int i, unsigned int j, unsigned int k) { unsigned int b; s1=i; s2=j; s3=k; b = ((s1 « 9) - s1) « 4; s1 = (s1 «4) - (b » 28); b = ((s2 « 2) - s2) « 3; s2 = (s2 « 3) - (b » 29); b = ((s3 « 6) - s3) « 1; s3 = (s3 « 1) - (b » 31); } double combT() { unsigned int b; b = ((s1 O, rozkładu wykładniczego. Jest to rzeczywiście łatw y rozkład, bo jeżeli U m a rozkład rów nom ierny U(0,1), to - ln U m a rozkład o gęstości g (patrz p. 3.1.1). Otrzym ujem y stałą c = V 2 e / n

i algorytm przyjm uje postać

A LG ORY TM 3.2 Repeat G eneruj X o rozkładzie wykładniczym E(0,1) G eneruj U o rozkładzie równomiernym U(0,1) until

Return X

W arunek zakończenia pętli pow tórzeń je st oczywiście rów noważny z w arunkiem (X - l)2 < - 2 ln U

W celu zakończenia generow ania zmiennej losowej o rozkładzie N (0,1) wystarczy teraz w yprodukować kolejne U o rozkładzie U(0,1) i zam ienić X na - X, gdy np. U < 0.5. M ożna ewentualnie skorzystać z generatora liczb losowych V o rozkładzie rów nom iernym U (-1,1) i zakończyć obliczenia działaniem X = XsignV. D alsze uproszczenie, polegające na tym, że generuje się tylko dwie zm ienne losow e V i X zam iast trzech U, V i X, opiera się na następującym lemacie: LEM A T 3.1. Jeżeli zm ienna losow a W m a rozkład równom ierny U(-1,1), to zm ienna losow a — IWI ma rozkład rów nom ierny U(0,1), a dw upunktow a zm ienna losow a sign W m a rozkład rów nom ierny na zbiorze {-1,1} i te dwie zm ienne losow e są niezależne.

Łatwy dowód tego lematu, ja k rów nież odpowiednią m odyfikację algorytmu, pozostawiam y Czytelnikowi. •

P rz y k ła d zasto so w an ia: ogon ro z k ła d u n o rm aln eg o Jako przykład praktycznego zastosow ania m etody eliminacji, a zarazem przykład różnych trików stosowanych w program ow aniu generatorów liczb losowych o zadanych rozkładach prawdopodobieństwa, om ówim y zagadnienie generow ania zmiennej losowej z prawego ogona rozkładu norm alnego N(0, 1), tzn. zmiennej losowej X o rozkładzie z gęstością f t (x) = ^ / ~ ~ * ^ ’ 2(1 - 0 ( 0 ) ’

x- t

(3 '7)

Przedstaw im y trzy algorytm y generow ania tej zmiennej. P rz y k ła d 2 (ogon rozkładu norm alnego I). G enerator liczb losowych X o rozkładzie z gęstością (3.7) m ożna skonstruować m etodą eliminacji, korzystając z oczywistej nierówności

x- t

e - 2/2 < t

Praw ą stronę tej nierówności m ożna przedstawić w postaci c(t) * g t(x), gdzie c(t) je st stałą (zależną od param etru obcięcia t) oraz g t(x) je st gęstością pewnego rozkładu prawdopodobieństwa. Żeby z funkcji x exp (-x2/2) zrobić gęstość praw dopodobieństw a na przedziale (t, ++OT)(x)

gdzie 1A(x), ja k zwykle, oznacza funkcję w skaźnikow ą zbioru A, tzn. funkcję przyjm ującą wartość 1 na zbiorze A oraz 0 poza nim. O graniczenie z góry gęstości (3.7) przyjm uje postać

[2 e~x2/2 f t (x) = J ------------------- < c (t) g t (x) V n 2(1 —0 ( t )) W 6 tW gdzie: c (t) = _ V ^ t (1 —® (t)) W standardowym algorytmie, konstruow anym m etodą eliminacji, generujem y X o rozkładzie z gęstością g t(x) oraz U o rozkładzie rów nom iernym U (0, 1) i sprawdzamy w arunek akceptacji c(t) Ugt(x) < f t (x) który po prostych przekształceniach przyjm uje postać UX < t

Jeżeli X m a rozkład praw dopodobieństw a o gęstości g t (x), to P {X < * } = I* g t (v)dv = 1 - exp ^ 1 ( t 2 - * 2)

w ięc zm ienną losow ą X m ożem y generow ać m etodą odw racania dystrybuanty: X = Vt 2 - 2 ln V , gdzie V je st zm ienną losow ą o rozkładzie rów nom iernym U(0, 1). Ta operacja w ym aga obliczania logarytm u naturalnego i pierw iastka kw adratow ego w głównej pętli algorytmu. Obliczanie logarytm u m ożna ew entualnie zastąpić generatorem zmiennej losowej o rozkładzie wykładniczym, natom iast obliczanie pierw iastka da się w yprow adzić poza głów ną pętlę algorytmu, gdzie będzie on obliczany tylko jeden raz dla każdej generowanej liczby losowej zam iast wielokrotnie, przy okazji każdego kandydata na tę liczbę. M ożna to wykonać w następujący sposób. W arunek { U X < t} je st tutaj rów noważny z warunkiem { U2X 2 < t2}. Ale X 2 = t2 -2lnV = 2(t2/2 - lnV). W prow adzając stałą r = t2/2 i zm ienną losow ą Y o rozkładzie w ykładniczym E(r, 1), otrzym ujem y algorytm A LG ORY TM 3.3 Repeat Generuj Y o rozkładzie wykładniczym E(r, 1) Generuj U o rozkładzie równomiernym U(0, 1) until 2U2Y < r Return x = 4 2 Y

Przykład 3 (ogon rozkładu normalnego II). W eźm y pod uw agę następującą nierówność: exp (-x2/2) < exp (t2/2 -tx),

x> t

Gęstość dom inująca w yraża się wzorem exp ( t2 / 2 g t (x) = — ------------- — — )— = t exp (t - tx)1(t,+»)(x) I" exp(t2 / 2 - ty)dy Algorytm przyjm uje postać: generuj X o rozkładzie z gęstością g t(x) oraz U o rozkładzie równom iernym U (0, 1) dopóty, dopóki nie zostanie spełniony w arunek U exp (t2 /2 - tx) < exp (-x2/2) czyli w arunek (X - t)2 < -2ln U

Pozostaje w ybór m etody generow ania zmiennej losowej X o rozkładzie z gęstością g t (x).

LEM A T 3.2. Jeżeli zm ienna losow a Y m a rozkład w ykładniczy E(0,1), to zm ienna losow a = t + Y/t ma rozkład o gęstości gt(x).

X

Po podstaw ieniu w w arunku akceptacji Y/t zam iast X - t oraz zmiennej losowej W o rozkładzie w ykładniczym E (0,1) zam iast - ln U, otrzym ujem y

A LG ORY TM 3.4 Repeat Generuj Y o rozkładzie wykładniczym E(0,1) Generuj W o rozkładzie wykładniczym E(0,1) until Y2 < 2t2W Return X = t + Y/t

W przykładach 2 i 3 mam y stałą c(t) = ^(t) / (t(l - O(t))), czyli efektywność eff(t) = l /c(t) m ierzona praw dopodobieństw em akceptacji, je st rów na odpowiednio: eff(l) = 0.6556, eff(2) = 0.8427, eff(3) = 0.9138, eff(4) = = 0.9470, itd. Algorytm przedstawiony w następnym przykładzie je st bardziej efektywny, ale w ym aga nieco dłuższego czasu przygotowawczego, potrzebnego do obliczenia stałych algorytmu. P rz y k ła d 4 (ogon rozkładu norm alnego III). Tak samo ja k w przykładach 2 i 3, generujemy zm ienną losow ą X o rozkładzie z gęstością (3.7), ale za gęstość dom inującą przyjm ujem y teraz gęstość przesuniętego rozkładu w ykładniczego g t—(x) = Ae~—x - \

x >t

D ystrybuantą tego rozkładu je st G a (x) = —- exp(-A ( x - 1) ) . G enerowanie zmiennej losowej X o takim rozkładzie m ożna wykonać m etodą odw racania dystrybuanty za pom ocą w zoru X = t - (1/ A)ln(1 - U ) gdzie U je st zm ienną losow ą o rozkładzie rów nom iernym U (0,1) lub w zoru X = t - Y / A , gdzie Y 1 je st zm ienną losow ą o rozkładzie w ykładniczym E(0,1). Zgodnie z ogólną zasadą metody eliminacji będziem y teraz generow ać niezależne zm ienne losow e U i V o rozkładzie równom iernym U(0,1) i obliczać X = t - ( l / A )ln U dopóty, dopóki nie zostanie spełniony w arunek cV g— ( X ) < f (X ) przy odpowiednio wybranej stałej c. D la ustalonych wartości param etrów t oraz X optym alna stała c = c(t, A) wynosi

M aksim um po prawej

stronie tego w zoru je st osiągane dla takiej wartości x > t, która

m aksym alizuje w ielkość A (x - 1) 2 - 1 x 2, czyli dla x = X, gdy X > t lub dla x = t w przeciwnym przypadku. Otrzym ujem y stąd optym alną w artość stałej c(X, t), w yrażającą się wzorem f ;2 exp — - A t 2 2A(1 - Q (t)) V 27 n

c(A, t ) =
t A< t

2A(1 - 0 ( t )) Param etr Xje st do naszej dyspozycji, m ożem y zatem przypisać mu taką wartość, żeby przy danym t zm inim alizować c(X, t). Tą w artością je st A opt = A ( t ) = —

(/

i stąd otrzym ujem y optym alną w artość copt=c(t):

a 12 + 4

+t)

(3.8)

V 27 n c (t) =

exp

( 4 12 + 4 - 3t)(Vt 2 + 4 + 1)

( 4 t 2 + 4 + 1)(1 - 0 ( t ))

8

(3.9)

Po podstaw ieniu w w arunku akceptacji optym alnych wartości param etrów X i c, po prostych przekształceniach w arunek ten przyjm uje postać V < exp ^ - 1 ( X - A ) 2 skąd ostatecznie otrzym ujem y algorytm A LG ORY TM 3.5 Repeat Generuj Y} o rozkładzie wykładniczym E(0, 1) Generuj Y2 o rozkładzie wykładniczym E(0, 1) X =t+Y /A until 1 - ( X - A ) 2 < Y2 Return X

Efektyw ność eff(t) = l/c(t), na mocy w zoru (3.9), wynosi eff(l) = 0.8765, eff (2) = 0.9336, eff(3) = 0.9609, eff(4) = 0.9749, itd. .

M eto d a elim in acji d la gęstości f (x) = cp(x)q(x) Podam y kilka w ersji m etody eliminacji, związanych ze szczególnymi postaciami gęstości. 1. Jeżeli gęstośćf m ożna przedstawić w postaci f (x) = cp(x)q(x), gdziep(x) je st także gęstością pewnego rozkładu, a funkcja q(x) przyjm uje wartości tylko w przedziale [0,1], to algorytm przybiera postać

A LG ORY TM 3.6 Repeat Generuj X o rozkładzie z gęstością p Generuj U o rozkładzie równomiernym U(0,1) until U < q(X) Return X 2.Jeżeli gęstość f m ożna przedstawić w postaci f (x) = cp(x)q(x), gdzie p(x) oraz q(x) są, odpowiednio, gęstością i dystrybuantą pew nych rozkładów prawdopodobieństwa, to algorytm przyjm uje postać

A LG ORY TM 3.7 Repeat Generuj X o rozkładzie z gęstością p Generuj Y o rozkładzie z dystrybuantą q until Y X Return X 5.Jeżeli gęstośćf m ożna przedstawić w postaci f (x) = cp(x) (l - q(t(x))), g d ziep(x) oraz q(x) są, odpowiednio, gęstością i dystrybuantą pewnych rozkładów prawdopodobieństwa, to algorytm przyjm uje postać

A LG ORY TM 3.10 Repeat Repeat G eneruj X o rozkładzie z gęstością p G eneruj Y o rozkładzie z dystrybuantą q until Y > t(X) Return X

P rz y k ła d 5 (dodatnia połów ka rozkładu norm alnego N(0,1)). D la dodatniej połówki rozkładu norm alnego N (0,1) mamy f (x ) = V2 / n e x p (-x 2 /2 ) , x > O, co m ożna zapisać w postaci

f (x) = j 2 e e n

(l —(1 —e -0

f (x) - a _ e fVni Przedstaw iam y tę gęstość w postaci

(3.18)

f ( x) = «1f 1( i ) g l (x) + «2f 2(x) g 2( x) gdzie 1 «2 = fl f i (x) = ■ 0

dla

0x ’,

0 < x < 1,

ci - 0

przy czy m M < ro . W tedy ^ M c /(i +1) = 1 i otrzym ujem y prosty algorytm: 1. 2.

W ygenerować indeks I e {1,2,...} w edług rozkładu prawdopodobieństw a P { I = i} = c-J(i + 1). D la danej w artości I = i w ygenerow ać zm ienną losow ą o rozkładzie z gęstością (i + l)x;, 0 < x < 1.

2 Zm ienną losow ą o rozkładzie z gęstością (i + l)x , 0 < x < l, m ożna łatw o wygenerow ać m etodą odwracania dystrybuanty albo za pom ocą w zoru X=max{Ui,U2,...,Uj+_i} gdzie Ul, ^ , n i są niezależnym i zm iennymi losowymi o jednakow ym rozkładzie równomiernym U(0,1). M oże się okazać, że w „dobrym" w ielom ianie aproksymującym gęstość nie w szystkie współczynniki ci są dodatnie. Zapiszm y ci, w postaci ci = ci+ - ci", W tedy

ci+, ci" > 0

f(x)=s cixi -1/2 N = [X] Generuj U o rozkładzie równomiernym U(0, 1) until a —5 X + ln — ----------------------— < k + N ln A —ln(N !) (1 + e x p ( a - 5 X ) ) 2 R etu rn N Istnieje w iele innych m etod generow ania zmiennej losowej o rozkładzie Poissona. K ilka algorytm ów , o których tutaj nie m ów iliśm y, m ożna znaleźć w cytowanych wyżej dwóch pracach A tkinsona oraz w dużej monografii D evroye'a (1985). R o zk ład geom etryczny Z m ienna losow a X m a rozkład geom etryczny G(p), jeżeli P{X = x} = (1 - p ) p x,

x = 0,1,2,...

N ajprostszy algorytm generow ania zm iennej losowej o tym rozkładzie je s t oparty na następującym lem acie:

LEMAT 3.10. Jeżeli zm ienna losow a X m a ro zk ła d w ykładniczy o g ęsto ści a e-ax , to zm ienna losow a [X ] m a ro zk ła d g eo m etryczn y G(e~a). N a mocy tego lematu, w celu w ygenerow ania zmiennej losowej X o rozkładzie (3.28) wystarczy w ygenerować zm ienną losową U o rozkładzie rów nom iernym U(0, 1) i obliczyć X = [ln U / ln p].

3.2.2. Rozkład wykładniczy

G ęstość praw dopodobieństw a rozkładu w ykładniczego E (0,X ) w yraża się wzorem

f e A x ) = ~ exP [^-

|

0< x U n < U n +1 i numerujemy je kolejnym i liczbam i O,l,2,... N um er pierwszej serii, w której n (długość serii) jest liczbą nieparzystą, przyjm ujem y za część całkowitą generowanej liczby losowej X, natom iast w artość R 1 pierwszej liczby w tej serii - za jej część ułam kow ą. A by przekonać się, że generow ana w ten sposób liczba losow a X m a rozkład w ykładniczy £(0 ,1 ), rozpatrzm y następujące zdarzenia losowe: E n (u) = { u > U l > U 2 >■■■ > U n } A n (u ) = { u > U l > U 2 > . ..> U n < U n + l }

Z auw ażm y, że A n(u) = E n (u) - E n+1 (u) oraz że zdarzenie £ (1 ) polega na zaobserw ow aniu serii „nieparzystej". N iech B 0(u), B 1(u),... będzie ciągiem niezależnych zdarzeń losow ych, tak ich ja k zdarzenie B(u). W tedy P{m < X < m + u} = P {B 0c(1) n B ‘ (1) n ... n B m_1 (1) n B m {u)} gdzie B c jest zdarzeniem przeciwnym do zdarzenia B.

W w yniku obliczeń otrzym ujem y

P {E n (u)} =

i....... idu d u j > > 2> >nj i

u u

u

u n 2

... u

....du n

n!

un u n+1 P {A n(u)} = P {E n} - P {E n+i) = — - ­ n! (n +1)! ro ( u.j/2n+1 u 2n+2 ^ = 1 - e -u P {B (u )} = £ ( 2 n + 1 )! ( 2 n + 2 )!) n=0

P{m < X < m + u} = P { (B c(1))m}P{B(u)} = = (e -1) m(1 - e -u) = F (m + u) - F (u) gdzie F je s t dystrybuantą rozkładu w ykładn iczego (3.30). Generator liczb losow ych o rozkładzie wykładniczym je st jednym z najczęściej używanych generatorów - bądź bezpośrednio do produkcji liczb losowych o tym rozkładzie, bądź jako pewien element w generatorach liczb losow ych o innych rozkładach. W bogatej literaturze m ożna znaleźć kilkanaście innych generatorów liczb losow ych o rozkładzie w ykładniczym (patrz np. m onografia D evroye'a (1986)). Jako pew ną ciekawostkę odnotujmy, że n niezależnych zm iennych losow ych o jednakow ym rozkładzie w ykładniczym E(0, 1) m ożna w ygenerow ać za pom ocą (n-1) liczb losowych o rozkładzie równomiernym U(0, 1) i jednej liczby losowej Z o rozkładzie gam ma r(n , 1): jeżeli U1:n-1,...,U n-1:n-1 je st statystyką pozycyjną z rozkładu rów nom iernego U(0,1), to Z U 1:n-1,Z(Z2:n-1U1:n-1) ,... .,Z(1-Un-1:n-1) je st ciągiem niezależnych zm iennych losow ych o jednakow ym rozkładzie w ykładniczym E(0, 1).

3.2.3. Rozkład normalny

Gęstość praw dopodobieństw a rozkładu norm alnego N (p,a) w yraża się wzorem f (X) = - ^ eXP a d 2n

--

r

\2 \ ' a

)

Jeżeli zm ienna losow a X m a rozkład norm alny N(jU,o), to (X- r )/ o m a rozkład norm alny N (0,1); dalej będziem y zajmowali się w ięc tylko zm iennymi losow ym i o rozkładzie norm alnym N(0, 1). G ęstość praw dopodobieństw a tego ro zk ład u oznaczam y przez ę, a d y strybuantę przez F.

Z e w zg lęd u n a liczne i b ard zo ró żn o ro d n e zasto so w an ia ro zk ła d u opracow ano w iele algo ry tm ó w g en ero w an ia liczb lo so w y ch o ty m rozkładzie.

n o rm aln eg o

M eto d ę od w ra ca n ia d y stry b u a n ty d la ro z k ła d u n o rm a ln e g o p rz e d sta w iliśm y w p.

3.1.1. Je st to m e to d a d o k ła d n a i p rzy dobrej im p lem en tacji k o m p u tero w ej (sch em at H o m era!) n a ty le szybka, że w arto je j się p rzy jrzeć n a sw oim kom puterze. P o d staw o w a jej zaleta p o leg a n a tym , że u ży w a ty lk o jed n ej liczb y losow ej o ro zk ład zie ró w n o m iern y m U (0,1) (por. p. 3.1.1).

M etodę elim inacji z g ęsto ścią ro zk ła d u w y k ład n iczeg o ja k o g ęsto ścią d o m in u jącą p o d aliśm y w p rzy k ła d zie l w p. 3.1.2. In n y a lg o ry tm o p raco w an y w ed łu g tej m etody, oparty n a pew nej szczególnej fak to ry zacji gęstości, p o d aliśm y w ty m sam y m p u n k cie w p rzy k ła d zie 5. M etod ę su p erp ozycji rozk ładów ilu stro w a liś m y w p. 3.1.3 g e n e ra to re m liczb lo so w y c h o ro z k ła d z ie n o rm aln y m w p rz y k ła d a c h l i 2. M etod ę R O U d la ro z k ła d u n o rm a ln e g o p rz e d s ta w iliś m y w p. 3 .1 .4 w p rzy k ład zie 2. W aru n e k ak cep tacji w p o d an y m ta m alg o ry tm ie 3.18 m iał postać {U2 < e~X /2}.W punkcie 3.1.2 m ów iliśm y o w aru n k ach szybkiej ak ce p tacji i szybkiej elim in acji, o p arty c h n a od p o w ied n io d o b ran y ch n ie rów nościach. R o zw ażany teraz w aru n ek akceptacji m ożna zapisać w p ostaci { X < - 4 1n U }. K o rz y sta ją c np. z n ie ró w n o śc i

3 „ 1 2 — 2 u H— u < —ln u 2

2

r 1I ---u 0

która je st praw dziw a dla każdego dodatniego y. W oryginalnej propozycji Chenga (1977) jest Y =912. O statecznie otrzym ujem y następujący algorytm:

A LG ORY TM (3.35) ^ d =

1 . =■, V 2 a —1

1 c = a H— , d

j i/i b = a —ln 4

Repeat Generuj U o rozkładzie równomiernym U(0, 1) Generuj U2 o rozkładzie równomiernym U(0, 1) V=d •ln - Ł 1 —U X = a e v, Z = U 2U ,, R = b + cV —X i f a k c e p tu j = fa ls e th e n a k c e p tu j = (R > ln Z) u n til a kc ep tu j

Return X Jeżeli zm ienne losowe X 1,..., X n są niezależne i zm ienna losow a X m a ro zk ład r ( a i), to zm ienna losow a X =

X

m a ro zk ład gam m a r ( a ) , gdzie a = ^ ” a t Ten fakt może być

wykorzystany w następujący sposób: je ż e li a > l, to zm ien n ą lo so w ą o tym ro z k ła d z ie m ożna w y g en ero w ać jako sumę [a] niezależnych zmiennych losowych o rozkładzie wykładniczym E(0, 1) i jednej zm iennej losowej o rozkładzie r ( { a } ) .

3.2.5. Rozkład beta G ęstość praw dopodobieństw a rozkładu beta B(a,b) w yraża się w zorem 1

f^ b (x) = ^ T x a— 1(1 —x )b— 1, B (a , b)

0 O są param etram i rozkładu. Zw racam y uw agę na to, że zarów no rozkład B(a, b), stalą B(a, b), ja k i czasami samą zm ienną losow ą o tym rozkładzie, oznaczamy takim samym symbolem; z kontekstu zawsze w iadom o o co chodzi, a unikam y w ten sposób nadm iaru oznaczeń.

Jeżeli zm ien n a lo so w a X m a ro zk ład b eta B (a,b), to zm ien n a lo so w a l - X m a rozkład b eta B(b, a), w ięc b ędziem y rozw ażali ty lk o przypadek, gdy a < b. Trzy najprostsze algorytm y generow ania zm iennych losow ych o rozkładach b eta są oparte na trzech następujących lem atach:

LEMAT 3.11 (Przypadek całkowitych a oraz b). Statystyka pozycyjna Uk:n z n-elementowej próby z rozkładu równomiernego U(0,1) m a rozkład beta B (k,n —k + 1) LEMAT 3.12. Jeżeli r(a ) oraz r (b) są dwiema niezależnymi zmiennymi losowymi o rozkładach gam m a z param etram i odpowiednio równymi a i b, to zm ienna losow a r (a )/( r (a) + r (b)) m a rozkład beta B (a, b).

LEMAT 3.13 (Algorytm Jdhnka (1964)). Jeżeli U, V są niezależnym i zm iennym i losowym i o jednakow ym rozkładzie równomiernym U(0,1) oraz a, b są dodatnimi stałymi, to rozkład warunkowy zm iennej losowej U 1a

U 1/a

+

V 1b

przy warunku {U1/a + V1/b < 1}, je s t rozkładem beta B(a,b). Ze w zględu na czas generowania liczb losowych, pierwszy z tych algorytm ów może konkurow ać z dw om a pozostałym i tylko w przypadku m ałych całkowitych w artości a oraz b, a drugi tylko wtedy, gdy dysponujem y szybkim i generatoram i liczb losow ych o rozkładach gam m a. A lgorytm Jdhnka je st efektywny w przypadku m ałych wartości param etrów a i b: r r praw dopodobieństw o spełnienia podanego w nim w arunku je s t rów ne (a+1) (b+1) T(a+b+1) W literaturze m ożna znaleźć kilkanaście szybkich i łatw ych do im plem entacji algorytm ów rozkładu beta. Ze w zględu na m etodykę konstrukcji, w arto zw rócić uw agę na algorytm zaproponowany przez Chenga (1978), oparty na następującym lemacie:

LEMAT 3.14 Jeżeli zm ienna losow a Y m a ro zkła d beta H rodzaju o gęstości y a-1 0

(3.41)

to zm ien n a lo sow a X — 7/(1 + Y) m a r o zk ła d beta B(a, b).

Zm ienną losow ą Y o rozkładzie (3.41) m ożna generow ać m etodą eliminacji z gęstością dom inującą 1 -1 ( \ 1dyA g *J y ) = E T + T 7 ’

r\ y >0

O dpowiednia dystrybuanta ma postać y1 G 2 u (y ) = —- — r , ^ (y) p + y 1 Z kolei

y >0 y

zm ienną losow ą Y o takiej dystrybuancie łatw o otrzym uje się m etodą odw racania

dystrybuanty. Standardow y algorytm przyjm uje postać: generow ać Y o rozkładzie z dystrybuanta Gx,p(y) i zm ienną losow ą U o rozkładzie rów nom iernym U (0, 1) dopóty, dopóki nie zostanie spełniony w arunek { c U g ^ (Y ) < 9 ab(Y)} i obliczyć X = Y /(1+Y ). Stala c pow inna oczyw iście zależeć od X oraz p,, a te param etry chcielibyśm y w ybrać tak, żeby ta stalą była ja k najmniejsza. W ybór takich optym alnych X i p zależy od param etrów a i b rozkładu

generowanej zmiennej losowej, ale - ja k zwykle - w olim y tu zdecydować się na jakieś proste, uniw ersalne (dla dużego zakresu w artości param etrów a i b) i łatw e do zaprogram ow ania procedury w yznaczania X oraz p, niż na rozw iązyw anie odpow iedniego zadania optym alizacji przy każdym odw ołaniu się do generatora. O ryginalną propozycją C henga je st

a X = \ 2ab —(a + b) 11 a + b —2 oraz

gdy a < 1 gdy a > 1

fa u = I— * lb

(Przypom inam y że zgodnie z przyjętymi na początku p.3.2.5 oznaczeniami, a < b.) W tedy 4 a abb . i \ a+b XB(a, b)(a + b )a

^

(3-42)

P o w y k o n a n iu o d p o w ie d n ic h p rz e k s z ta łc e ń o trz y m u je m y a lg o ry tm ALGORYTM 3.36 a =a +b 12ab —a I f a < 1 then R =a else R = J ----------V a —2 Y= a +R Repeat Generuj U o rozkładzie równomiernym U(0, 1)

V =- l n - Ł W = aeV P 1 —U Generuj U2 o rozkładzie równomiernym U(0, 1) a until a ln --------- + vV —ln 4 > ln(U,2U? ) b+W Return X = ■ b+Y Stalą (3.42.) nigdy nie przekracza liczby 4, a jeżeli a > l, to ta stała je st nie w iększa od 4le ~ l .47.

3.2.6. Rozkład Cauchy’ego Gęstość praw dopodobieństw a rozkładu Cauchy'ego C(Q,X) w yraża się wzorem r r \ X 1 f ep (x ) = 7T X + (X —0)

—“ < x < + “

Jeżeli zm ienna losow a Y m a rozkład C auchy'ego C(0,X) to zm ienna losow a X = (Y - 0)/X ma rozkład Cauchy'ego C(0, 1) o gęstości f (x) =

— z, T 1+ X

—ro < x < ( 3 . 4 3 )

w ięc w dalszym ciągu będziem y zajmowali się tylko rozkładem C(0,1) o gęstości (3.43). D ystrybuanta rozkładu C(0.1) m a postać

F (x) = — arctg (x) + 1 n 2

(3.44)

skąd w ynika prosty sposób generow ania liczb losowych o tym rozkładzie m etodą odwracania dystrybuanty; X = tg(Un), gdzie U je st zm ienną losow ą o rozkładzie rów nom iernym U

1 , 1 jj .

Przedstaw im y szybki i elegancki algorytm generow ania liczb losow ych o rozkładzie Cauchy'ego C (0 ,1) za pom ocą pewnej kom binacji m etody superpozycji rozkładów i m etody elim inacji. LEMAT 3.15. Jeżeli zm ien n a losow a X m a ucięty ro zk ła d C auchy'ego C (0,1) o g ęsto ści 2

f u (x) = i [0 ,

1

n 1+ x2 p o za tym

g dy

,

- 1< x < 1

(3 45)

to zm ien n a lo sow a Y, któ ra z p ra w d o p o d o b ień stw em 1/2 j e s t id en tyczn a ze zm ienną losową X oraz z praw dopodobieństw em 1 /2 je s t identyczna ze zm ien n ą losow ą 1/X, m a ro zk ła d C auchy'ego C(0,1). pb 2 Dowód: Pam iętając o tym , że I d t /(1 + 1 ) = arctg (b) - arctg (a), dla zmiennej losowej Y i dla Ja wartości argum entu y < - 1 łatw o otrzymujemy:

P {Y < y} = 1 P{ X < y} + 1 p j - L < y U

i

0 + 1 P j 1 < X < o l = 1 - p0 dt 2 : t 2 [y J 2 n Jl/y j1/y -11 = 1 1 = — arctgy + n 2 czyli dystrybuantę (3.44). Podobnie dowodzim y lematu, gdy y € (-1,1) oraz gdy y >

2A

i 1 1 1 -1

f 'max

-----------’ ' ------] 'Fftin ¡in o * ? I I I I 0 1 x -

Rys. 3.8 Zm ienną losow ą X o rozkładzie z gęstością f u(x) w ygenerujem y m etodą AC (patrz algorytm 3.15), przyjm ując za gęstość dom inującą g gęstość rozkładu rów nom iernego U(-1, 1) oraz z a f gęstość stałą (używamy tu oznaczeń ja k w algorytm ie 3.15). Poniew aż g(x) = Ą, w ięc zaf 1 przyjm ujem y f u(x) - f a x - 1/2) , gdzie f max = max f u (x). W tedy f 2 (x) = f max - Ą (patrz rys. 3.8) Otrzym ujem y następujący algorytm: A LG ORY TM 3.37

Generuj X o rozkładzie równomiernym U(-1, 1) Generuj U o rozkładzie równomiernym U(-1, 1) If

1 u >/

(x

-1 1

then Generuj X o rozkładzie równomiernym U(-1, 1) Return X Poniew aż f rmix = 2/n, w ięc dodatkowego generow ania X w edług rozkładu równom iernego U(-1, 1) przyjm uje postać (U + 0.27324)(1 + X 2) > 1.27324. ten w arunek je st spełniony z praw dopodobieństw em j " / (x)d x = 2 ( / mx - 1 / 2 ) = 0.27324

3.2.7. Rozkłady a-stabilne

R ozkłady a -stabilne tw orzą obszerną klasę, która zaw iera m.in. rozkład normalny (a = 2) i rozkład Cauchy'ego (a = 1). Kłopoty z generowaniem zm iennych losow ych o takich rozkładach polegają przede w szystkim na tym, że poza w yjątkow ym i sytuacjami nie je st znany w zór dla gęstości lub dystrybuanty tego rozkładu, a w iec takie m etody, ja k m etoda odw racania d y stry b u an ty lub m etoda elim inacji, n ie zn ajd u ją tutaj zastosow ania. Z m ienna losow a X m a ro zk ła d a-stabilny, 0 < a < 2. z param etrem skali o > 0, z param etrem asym etrii yff€[-1, 1] i z param etrem położenia p e R , jeżeli jej funkcja charakterystyczna ę(t) je st określona wzorem - a a |t| a f 1 - ifisignttg ~ ~ j + i^ t,

gdy

a * 1

ln ę ( t ) = - cr|t 1 1 + iPsignt — ln|t| J + i^ t

gdy

a =1

Ten rozkład oznaczam y przez Sa (o,P,p). Jeżeli zm ienna losow a X m a rozkład Sa (1,P,0), to zm ienna losow a Y określona w zorem oX +u, gdy Y =1 2 \a X H— ^ < r ln ^ + u , n

a* 1 gdy

a =1

m a rozkład Sa (o,P,p), w ystarczy w ięc skonstruow anie algorytm u generow ania zm iennej losow ej X o rozkładzie a-stab iln y m X a(1,P, 0). M ożna to zrobić w następujący sposób, podany po raz pierwszy w pracy Chambersa, M allow sa i Stucka (1976): wygenerow ać zm ienną losow ą U o rozkładzie równomiernym U(-n/2, n/2) oraz niezależnie zm ienną losow ą W o rozkładzie w ykładniczym £(0,1), a następnie: 1) jeżeli a = l, obliczyć r

A n y - W cos U n * P V W U - ^ ln -2- ---------+p u 2

1) jeżeli a Y l, obliczyć s in a ( u + B a p) f cos (U - a ( u + B a j y X = Sa.P •

(cos U )1/a

W

(1-a )/a

gdzie f B aj = a ^a rctg{^tg

f x l/(2a) S ap = ^ + p 2tg 2 ^

\

N ie podajem y tu żadnych dowodów; zwięzły w ykład na ten tem at oraz bardziej specjalistyczną bibliografię m ożna znaleźć w książce Janickiego i W erona (1994) oraz w pracach W erona (1996a. 1996b).

3.3. Związki między rozkładami Podstaw ą konstrukcji generatorów liczb losow ych o zadanych rozkładach praw dopodobieństw a są związki teoretyczne m iędzy różnymi zm iennymi losowymi. W bieżącym podrozdziale zam ieszczam y listę takich związków. Odnalezienie inform acji na tem at konkretnych, interesujących Czytelnika rozkładów, ułatwi schemat graficzny przedstawiony na rys. 3.9 1. R ozk ład n o rm aln y i ro zk ła d C a u c h y ’ego. Jeżeli zm ienne X i Y są niezależnymi zmiennymi losowymi o rozkładzie normalnym N(0,1), to zm ienna losow a Z = X /Y ma rozkład Cauchy'ego C (0, 1). 2. R ozk ład n o rm aln y i ro zk ła d L ap lace'a. Jeżeli X l,X2,X3,X4 są niezależnymi zmiennymi losowymi o rozkładzie normalnym N (0,1), to zm ienna X lX 4 - X 2X 3 m a rozkład Laplace'a z gęstością f (x) = 1 e x p ^—1 |x —#|j. 3. R ozk ład n o rm aln y i ro zk ła d ch i-k w ad rat. Jeżeli X i,....X n są niezależnymi zmiennymi losowymi o jednakow ym rozkładzie norm alnym N(0,1), to zm ienna losow a

X l ma

rozkład chi-kw adrat X 2 o n stopniach swobody. 4. R ozk ład n o rm aln y i ro zk ła d lo g n o rm aln y . Jeżeli zm ienna losow a X m a rozkład norm alny N (p,o), to zm ienna losow a Y = ex m a rozkład lognorm alny z param etram i (p,o).

Rys. 3.9

5.

R o zk ład n o rm aln y i ro z k ła d ró w n o m iern y .

(a) Jeżeli X i 7 są niezależnymi zm iennymi losowymi o rozkładzie rów nom iernym na przedziale (0,1), to U = V - 2 ln X • cos(2nY), V = V - 2 ln X • sin(2nY) są niezależnymi zmiennymi losowymi o rozkładach norm alnych N(0,1). (b) Jeżeli X i 7 są niezależnymi zm iennymi losowymi o rozkładzie norm alnym N(0,1), to zm ienne losow e U = X 2 + Y 2,

V = arc sin (x / VX 2 + Y 2 j

są niezależne, U m a rozkład

w ykładniczy £ (0,2), V m a rozkład rów nom ierny na przedziale (-n/2, n/2). 6. R o zk ład n o rm aln y , ro z k ła d c h i-k w a d ra t i ro zk ła d t-S tu d en ta. Jeżeli zm ienna losow a X m a rozkład norm alny N (0,1), zm ienna losow a 7 m a rozkład chikw adrat o n stopniach swobody i te zm ienne są niezależne, to zm ienna losow a

t = x / T Y T n ma rozkład Studenta (rozkład t-Studenta) o n stopniach swobody. 7. R ozk ład ró w n o m iern y i ro zk ła d C au ch y 'eg o . Jeżeli zm ienna losow a X m a rozkład równom ierny na przedziale (-n/2, n/2), to zm ienna losow a 7 = tg X m a rozkład Cauchy'ego C(0,1). 8 . R ozk ład ró w n o m iern y i ro zk ła d P areto . Jeżeli zm ienna losow a U m a rozkład równom ierny U(0,1), to zzm ienna losow a X = b U 1a m a rozkład Pareto o dystrybuancie

Fab (x) = 1 - (b / x )a, x > b. 9. R ozk ład ró w n o m iern y i ro zk ła d beta. (a) Jeżeli X 1v.., X n są niezależnymi zmiennymi losowymi o rozkładzie równom iernym na przedziale (O, 1), to k-ta statystyka pozycyjna X k:n m a rozkład b eta B (k,n - k + 1), a zm ienna losow a X l:n - X k:n (dla l > k) m a rozkład beta B (l - k,n - l+ k +1).

(b) Jeżeli C ,V są niezależnymi zm iennymi losowymi o jednakow ym rozkładzie równomiernym U(0, 1) oraz a, b są dodatnimi stałymi, to rozkład w arunkowy zmiennej losowej U 1a U 1/ a + V 17b przy w arunku { U l/ a + V 1/b < 1}, je st rozkładem beta B(a,b). 10. R o zk ład ró w n o m iern y i ro zk ła d logistyczny. (a) Jeżeli zm ienna losow a X m a rozkład rów nom ierny na przedziale (0,1), to zm ienna losow a Y = a - bln (1/X - 1) m a rozkład logistyczny L(a, b) o gęstości exp ( - (x - a ) / b ) b(l + exp( - ( x - a) / b))2 ’

b>Q

(b) Jeżeli zm ienna losow a U m a rozkład rów nom ierny U(0,1), to zm ienna losowa X = ln (U /(1 - U)) m a rozkład logistyczny o dystrybuancie F (x) = 1/ (1+e’x). 11.R ozkład ró w n o m iern y i ro zk ła d potęgow y. (a) Jeżeli X i,...., X k są niezależnymi zm iennymi losowymi o równomiernym na przedziale (0,1), to max (X j,...., Xk) je st zm ienną losow ą o gęstości f(x ) = kxk_1, 0 < x 0, j = 1,2,... , m} a pow ierzchnie takich w ielościanów - przez przekształcenia zbioru

= {(xŁ,...,xm) : i m ^ = 1,xy > 0, j = 1,2,...,m } . Jest oczywiste, że jeżeli

X je st zm ienną losow ą o rozkładzie równomiernym na zbiorze Q oraz F je st nieosobliwym liniowym przekształceniem przestrzeni R m, to F (X .) je st zm ienną losow ą o rozkładzie równom iernym na obrazie F (Q ) zbioru Q w tym przekształceniu.

4.2.2. Rozkład równomierny na sferze i na kuli w R m Przypuśćmy, że zm ienna losow a X = (X 1, ^ , X m) m a rozkład rów nom ierny na sferze

Sm=| (

x

m): I x2 = i [ j =1

(« )

Zauważmy, że jeżeli A je st m acierzą ortonorm alną (tzn. m acierzą pewnego obrotu w R m) i jeżeli X ma rozkład równom ierny na sferze Sm, to zm ienna losow a AX m a rów nież rozkład rów nom ierny na tej sferze. M ówimy, że m -w ym iarowa zm ienna losow a Z = (Z1,...,Zm) m a rozkład sferycznie konturowany, jeżeli gęstość g z( z ^ . . ,,zm) jej rozkładu zależy tylko od llzll2 II II = Y m =_1 zj 2 Najprostszym przykładem takiego rozkładu je st łączny rozkład m niezależnych zm iennych losow ych o rozkładach norm alnych N(U,1):

gz (zl ,..., Zm) = ( 2 n ) - m/2 eXP | - 1 NI2 1 2

2

Poniew aż dla każdej ortonormalnej macierzy A mam y ||Az|| = ||z|| , w ięc jeżeli Z ma rozkład sferycznie konturowany, to A Z m a taki sam rozkład jak Z. Prowadzi to do następującego algorytmu generow ania zmiennej losowej X o rozkładzie równomiernym na sferze Sm:

ALG ORY TM 4.1 Wygeneruj m niezależnych zm iennych losowych Z i,..., Z m o rozkładzie normalnym N(0,1) O b liczX j = Zj/||z||, j = 1,2,..., m Return X = (Xi,...,Xm) Inna m etoda generow ania zmiennej losowej X oparta je st na następującym lemacie:

L e m a t 4.1. Jeżeli Z i,...,Z k m a rozkład równom ierny na k-wym iarowej sferze Sk , R je s t zm ienną losową o rozkładzie z gęstością .k-1 cr h (r ) =

1- r ^ 0,

je że że 0 < r ą j = l ,2 ,..., m } j =1 Z lem atów wynika, że generow anie punktów losow ych na m-wym iarowym sympleksie lub na jego powierzchni sprowadza się do sprawnego w yznaczenia statystyk pozycyjnych w m- lub (m l )-elem entowym ciągu niezależnych zm iennych losowych o rozkładzie równomiernym U(0,1). N ajprostszy algorytm generow ania ciągu U 1:m, . . ,,U m:m polega oczywiście na uporządkow aniu ciągu U i,...,Um w ciąg niemalejący, jednak znane są bardziej efektywne algorytmy. Oto dwa z nich: ALG ORY TM 4.5 Wygeneruj ciąg E 1, . , wykładniczym E(0,1) S =1 ^

;

E m+1 i niezależnych zm iennych losowych o rozkładzie

U„m = 0

F or j = l to m do Uj:m = Uj-1:m +Ej /S Return U1:m,-,Um:m ALG ORY TM 4.6 Um+1:m 1 F or j = m downto 1 do G eneruj U o rozkładzie równomiernym U(0,1) Uj:m = U UJ Uj+1:m Return U 1 :m,.,Um:m Algorytm 4.5 je st oparty na następującym lemacie: LEM A T 4.4. N iech U1:m,...,U m:m będzie ciągiem statystyk pozycyjnych ciągu U1,...,U m niezależnych zm iennych losowych o rozkładzie równomiernym i niech Uo:m = 0 oraz Um+1:m = 1. Zmienne losowe Uj:m - Uj-1:m, j= l, 2,..., m+1, m ają taki sam rozkład ja k zm ienne losowe E 1/S,E2/S,...,E m+1/S, gdzie E 1,... ,E m+1 je s t ciągiem niezależnych zmiennych losowych o rozkładzie wykładniczym E(0,1) oraz S = W ynika stąd prosty algorytm generow ania punktów o rozkładzie równom iernym na powierzchni sympleksu: wygeneruj ciąg E 1,...,Em niezależnych zmiennych losowych o rozkładzie wykładniczym E(0,1) i oblicz E 1 /S, E 2/S,..., Em/S, gdzie S = E 1 + E 2 + ... + Em. Algorytm 4.6 je st oparty na następującym lemacie: Lem at 4.5. N iech U1:m,...,U m:m będzie ciągiem statystyk pozycyjnych ciągu U1,...,U m niezależnych zm iennych losowych o rozkładzie równomiernym. D la każdego k = 0,l,...,m-l, zmienne

losowe

Um:m,Um1 :m,...,Uk:m,

mają

taki

sam

rozkład

ja k

zmienne

losowe

I/n t 7-1/nr 7-1/(m— 1) t 7-1/nj r1/(m— 1) t j\!(m—k) U n ,U n U n—1 ,•••, U n U n—1 "'U m— k 7- 7-

4.3. Wielowymiarowy rozkład normalny Prosta m etoda Boxa-M ullera (1958) generow ania dwuwymiarowej zmiennej losowej o rozkładzie normalnym je st oparta na następującym spostrzeżeniu: jeżeli Ui i U2 są dwiem a niezależnymi zm iennymi losowymi o jednakow ym rozkładzie równomiernym U(0,1), to zm ienne losowe X 1 i X 2 określone wzoram i X = yj—2 ln (U ) c o s (2 ^ 7 2)

X

= *J—2 ln (U ) sin ( 2 ^ U 2 )

są niezależne i m ają jed n a ko w y rozkład norm alny N(0,1). Z a pom ocą tej m etody otrzym ujem y jednocześnie dwie niezależne zm ienne losowe o jednakow ym rozkładzie normalnym N(0, 1). W zastosow aniach często trzeba korzystać z w ielow ym iarow ych rozkładów normalnych, w których poszczególne zm ienne losow e są ze sobą skorelowane. N iech X A =

CT1,2 •••CT1,n

CT2, CT2,2 •••CT2,n

_CTn, CTn,2 ••• n,n będzie m acierzą kowariancji n-wymiarowej zmiennej losowej X = (X 1,...,X n) o rozkładzie normalnym, takim że w artość oczekiw ana E X i = 0 dla każdego i = 1,2, ... ,n. Jeżeli Z = (Z1,..., Z n) oraz w szystkie Z , i= 1,2, ... , n, są niezależne i m ają taki sam rozkład norm alny N(0,1), to zm ienna losow a C Z , gdzie C je st pew ną m acierzą nieosobliwa, m a n-w ym iarow y rozkład norm alny z T T m acierzą kowariancji C C (C je st m acierzą transponow aną macierzy C). W celu w ygenerow ania n-wymiarowej zmiennej losowej X z daną m acierzą kowariancji A wystarczy w ięc skonstruować odpow iednią macierz C, taką żeby C C = A, wygenerow ać n niezależnych zm iennych losowych Z 1, ...,Zn o jednakow ym rozkładzie normalnym N (0, 1) i obliczyć X = CZ. M acierz C m ożna skonstruować w następujący sposób: i,1

i,1

,1 1/2 i—1 ci,i = X —X C2 r=1 ) V f J—1 ci, j = C1 X j —X c i,rc J,r V r=1

f

Ci, j = o,

g dy

gdyi > J

i