Ręczne generowanie liczb losowych

Jak mogę ręcznie wygenerować liczbę losową z danego rozkładu, na przykład 10 realizacji ze standardowego rozkładu normalnego?

Jeśli „ręcznie” obejmuje „mechaniczne”, masz wiele dostępnych opcji. Aby symulować zmienną Bernoulliego z połową prawdopodobieństwa, możemy rzucić monetą: dla ogonów, 1 dla głów. Aby zasymulować rozkład geometryczny, możemy policzyć, ile rzutów monetą potrzeba, zanim zdobędziemy głowy. Aby zasymulować rozkład dwumianowy, możemy podrzucić naszą monetę n razy (lub po prostu podrzucić n monet) i policzyć głowy. „Quincunx” lub „maszyna fasoli” lub „Galton skrzynka” jest bardziej kinetycznej alternatywę - dlaczego nie ustawić jeden do działania i zobaczyć na własne oczy ? Wygląda na to, że nie ma czegoś takiego jak „moneta ważona”$0$$1$$n$$n$ale jeśli chcemy zmienić parametr prawdopodobieństwa naszej zmiennej Bernoulliego lub zmiennej dwumianowej na wartości inne niż , igła Georges-Louis Leclerc, Comte de Buffon pozwoli nam to zrobić. Aby zasymulować dyskretny rozkład równomierny na { 1 , 2 , 3 , 4 , 5 , 6 } , rzucamy sześciościenną kostką. Fani gier fabularnych napotkają bardziej egzotyczne kości , na przykład kości czworościenne, aby próbkować jednolicie od { 1 , 2 , 3 , 4 $p = 0.5$$\{1, 2, 3, 4, 5, 6\}$$\{1,2,3,4\}$, podczas gdy za pomocą pokrętła lub ruletki można iść jeszcze dalej. ( Kredyt na zdjęcia )

Czy musielibyśmy być wściekli, aby generować losowe liczby w ten sposób dzisiaj, gdy jest tylko jedno polecenie na konsoli komputera - lub, jeśli mamy dostępną odpowiednią tabelę liczb losowych, jedną wyprawę do bardziej zakurzonych rogów półki na książki? Być może, choć w doświadczeniu fizycznym jest coś przyjemnego w dotyku. Ale dla osób pracujących przed erą komputerową, a właściwie przed szeroko dostępnymi wielkoskalowymi tabelami liczb losowych (których więcej później), ręczne symulowanie zmiennych losowych miało bardziej praktyczne znaczenie. Kiedy Buffon zbadał paradoks petersburski- słynna gra polegająca na rzucaniu monetami, w której suma wygranych gracza podwaja się za każdym razem, gdy rzucane są głowy, gracz przegrywa z pierwszymi ogonami, a spodziewana wypłata jest nieskończenie intuicyjna - musiał zasymulować rozkład geometryczny za pomocą . Aby to zrobić, wydaje się, że wynajął dziecko, aby rzuciło monetą, aby zasymulować 2048 odtworzeń gry w Petersburgu, rejestrując liczbę rzutów przed zakończeniem gry. Ten symulowany rozkład geometryczny jest odtworzony w Stigler (1991) :$p=0.5$

Tosses Frequency
1      1061
2      494
3      232
4      137
5      56
6      29
7      25
8      8
9      6

W tym samym eseju, w którym opublikował to empiryczne dochodzenie w sprawie paradoksu petersburskiego, Buffon przedstawił także słynną „ igłę Buffona ”. Jeśli płaszczyzna jest podzielona na paski równoległymi liniami w odległości od siebie, a igła o długości l ≤ d jest zrzucana na nią, prawdopodobieństwo, że igła przecina jedną z linii, wynosi 2 $d$$l \leq d$ .$\frac{2l}{\pi d}$

Igła Buffona może zatem służyć do symulacji zmiennej losowej lubX∼Dwumianowy(n,2$X \sim \text{Bernoulli}(\frac{2l}{\pi d})$, i możemy dostosować prawdopodobieństwo sukcesu, zmieniając długości naszych igieł lub (być może wygodniej) odległość, w której rządzimy liniami. Alternatywnym zastosowaniem igieł Buffona jest niesamowicie nieefektywny sposób na znalezienie probabilistycznego przybliżeniaπ. Obraz (kredyt) pokazuje 17 zapałek, z których 11 przecina linię. Gdy odległość między liniami liniowymi jest równa długości zapałki, tak jak tutaj, oczekiwany odsetek przekroczenia zapałek wynosi$X \sim \text{Binomial}(n,\frac{2l}{\pi d})$$\pi$ , a tym samym można oszacować π w podwójnej odwrotności obserwowanego frakcji: tu uzyskać π =2⋅$\frac{2}{\pi}$$\hat \pi$. W 1901 Mario Lazzarini osiągając przeprowadzili doświadczenia stosując 2,5 cm igieł z liniami 3 cm jedna od drugiej, i po otrzymano 3408 rzutów gatunku =$\hat \pi = 2 \cdot \frac{17}{11} \approx 3.1$ . Jest to dobrze znane racjonalne doπ, z dokładnością do sześciu miejsc po przecinku. Badger (1994) dostarcza przekonujących dowodów na to, że było to oszukańcze, zwłaszcza że aby być pewnym 95% pewności co do sześciu miejsc dziesiętnych za pomocą aparatu Lazzariniego, należy rzucić 134 biliony igieł porzucających cierpliwość! Z pewnością igła Buffona jest bardziej przydatna jako generator liczb losowych niż jako metoda szacowaniaπ.$\hat \pi = \frac{355}{113}$$\pi$$\pi$

Jak dotąd nasze generatory były rozczarowująco dyskretne. Co jeśli chcemy symulować rozkład normalny? Jedną z opcji jest uzyskanie losowych cyfr i użycie ich do utworzenia dobrych dyskretnych aproksymacji do równomiernego rozkładu na , a następnie wykonanie pewnych obliczeń w celu przekształcenia ich w losowe odchylenia normalne. Pokrętło lub koło ruletki może podawać cyfry dziesiętne od zera do dziewięciu; kostka może generować cyfry binarne; jeśli nasze umiejętności arytmetyczne poradzą sobie z zabawniejszą bazą, zrobiłby to nawet standardowy zestaw kości. Inne odpowiedzi obejmowały bardziej szczegółowo takie podejście oparte na transformacji; Odraczam dalszą dyskusję na ten temat do końca.$[0,1]$

Pod koniec XIX wieku użyteczność rozkładu normalnego była dobrze znana, dlatego też statystycy chętnie symulowali przypadkowe odchylenia normalne. Nie trzeba dodawać, że długie obliczenia ręczne nie byłyby odpowiednie, gdyby nie przede wszystkim skonfigurować proces symulacji. Po ustaleniu generowanie liczb losowych musiało być stosunkowo szybkie i łatwe. Stigler (1991) wymienia metody stosowane przez trzech statystyków tej epoki. Wszyscy badali techniki wygładzania: losowe odchylenia normalne były oczywistym przedmiotem zainteresowania, np. W celu symulacji błędu pomiaru, który musiał zostać wyrównany.

Wybitny amerykański statystyk Erastus Lyman De Forest był zainteresowany wygładzeniem tabel życia i napotkał problem, który wymagał symulacji wartości bezwzględnych odchyleń normalnych. W tym, co okaże się działającym motywem, De Forest naprawdę pobierał próbki z pół-normalnej dystrybucji . Co więcej, zamiast używać standardowego odchylenia jednego ( jesteśmy przyzwyczajeni do nazywania „standardowym”), De Forest chciał „prawdopodobnego błędu” (odchylenie środkowe) jednego. Taka była forma podana w tabeli „Prawdopodobieństwo błędów” w dodatkach do „Podręcznika astronomii sferycznej i praktycznej, tom II”$Z \sim N(0, 1^2)$William Chauvenet . Z tej tabeli De Forest interpolował kwantyle o rozkładzie półnormalnym, od do p = 0,995 , które uważał za „błędy o równej częstotliwości”.$p=0.005$$p=0.995$

Jeśli chcesz zasymulować rozkład normalny, po De Forest, możesz wydrukować tę tabelę i pokroić ją. De Forest (1876) napisał, że błędy „zostały wpisane na 100 kawałków kart o równej wielkości, które zostały wstrząśnięte w pudełku i wszystkie wyciągnięte jeden po drugim”.

Astronom i meteorolog Sir George Howard Darwin (syn przyrodnika Charlesa) zmienił sytuację, opracowując coś, co nazwał „ruletką” do generowania przypadkowych odchyleń normalnych. Darwin (1877) opisuje, jak:

$x$$\frac{720}{\sqrt \pi} \int_0^x e^{-x^2} dx$$+$$-$$+$$-$

„Indeks” należy tutaj czytać jako „wskaźnik” lub „wskaźnik” (por. „Palec wskazujący”). Stigler wskazuje, że Darwin, podobnie jak De Forest, używał pół normalnej skumulowanej dystrybucji wokół dysku. Kolejne użycie losowej monety w celu uzyskania pełnego rozkładu normalnego. Stigler zauważa, że ​​nie jest jasne, jak dokładnie została wyskalowana skala, ale zakłada, że ​​instrukcja ręcznego zatrzymania dysku w połowie obrotu była „w celu zmniejszenia potencjalnego odchylenia w kierunku jednej części dysku i przyspieszenia procedury”.

Sir Francis Galton , nawiasem mówiąc, pół kuzyn Karola Darwina, został już wspomniany w związku ze swoim quincunx. Podczas gdy mechanicznie symuluje to rozkład dwumianowy, który według twierdzenia De Moivre – Laplace'a jest uderzająco podobny do rozkładu normalnego (i jest czasami wykorzystywany jako pomoc dydaktyczna w tym temacie), Galton stworzył o wiele bardziej skomplikowany schemat, gdy chciał próbka z rozkładu normalnego. Jeszcze bardziej niezwykłe niż niekonwencjonalne przykłady na górze tej odpowiedzi, Galton opracował normalnie rozmieszczone kości- lub dokładniej, zestaw kostek, które dają doskonałe dyskretne przybliżenie do rozkładu normalnego z medianowym odchyleniem. Te kości, pochodzące z 1890 roku, są zachowane w kolekcji Galton w University College London.

W artykule z 1890 r. W Nature Galton napisał, że:

Jako narzędzie do losowego wybierania nie znalazłem nic lepszego niż kości. Najbardziej żmudne jest tasowanie kart pomiędzy kolejnymi losowaniami, a metoda mieszania i mieszania oznaczonych piłek w torbie jest jeszcze bardziej nużąca. Teetotum lub jakaś forma ruletki korzystne jest to, ale kości są lepsze niż wszystkie. Kiedy są wstrząśnięci i wrzucani do kosza, tak różnie walczą ze sobą i o żebra kosza, że ​​szaleją wokół, a ich pozycje na początku nie dają żadnej zauważalnej wskazówki, co będą po nawet pojedynczy dobry wstrząs i podrzucanie. Szanse, jakie daje kość, są bardziej zróżnicowane niż się powszechnie uważa; istnieją 24 równe możliwości, a nie tylko 6, ponieważ każda twarz ma cztery krawędzie, które można wykorzystać, jak pokażę.

$+$$-$$1 \frac 1 4$

Laboratorium eksperymentów statystycznych matematycznych Raazesh Sainudiin obejmuje projekt studencki z University of Canterbury, NZ, odtwarzający kości Galtona . Projekt obejmuje badanie empiryczne polegające na wielokrotnym rzucie kostką (w tym empiryczny CDF, który wygląda uspokajająco „normalnie”) oraz dostosowanie wyników kostek, aby były zgodne ze standardowym rozkładem normalnym. Wykorzystując oryginalne wyniki Galtona, istnieje również wykres dyskretnego rozkładu normalnego, za którym faktycznie podążają wyniki kości.

Na wielką skalę, jeśli jesteś przygotowany na rozciągnięcie „mechanicznego” na elektryczne, zauważ, że epicka RAND A Million Random Digits z 100 000 normalnych odchyleń była oparta na pewnego rodzaju elektronicznej symulacji koła ruletki. Z raportu technicznego (George'a W. Browna, pierwotnie z czerwca 1949 r.) Wynika:

W ten sposób zmotywowani ludzie RAND, przy pomocy personelu inżynieryjnego Douglas Aircraft Company, zaprojektowali koło do ruletki elektro-ruletkowej w oparciu o wariant zaproponowany przez Cecila Hastingsa. Na potrzeby tego wystąpienia wystarczy krótki opis. Źródło impulsu o losowej częstotliwości było bramkowane impulsem o stałej częstotliwości, mniej więcej raz na sekundę, zapewniając średnio około 100 000 impulsów w ciągu jednej sekundy. Obwody standaryzacji impulsów przekazywały impulsy na pięciopunktowy licznik binarny, tak że w zasadzie maszyna jest jak koło ruletki z 32 pozycjami, wykonując średnio około 3000 obrotów na każdy obrót. Zastosowano konwersję binarną na dziesiętną, wyrzucając 12 z 32 pozycji, a otrzymaną losową cyfrę wprowadzono do stempla IBM, uzyskując tabele kart perforowanych losowych cyfr.

$\chi^2$testy częstotliwości cyfr nieparzystych i parzystych wykazały, że niektóre partie wykazywały niewielki brak równowagi. W niektórych partiach było gorzej niż w innych, co sugeruje, że „maszyna przestała działać w ciągu miesiąca od momentu jej dostrojenia… Wskazania wskazują, że ta maszyna wymagała nadmiernej konserwacji, aby utrzymać ją w doskonałym stanie”. Znaleziono jednak statystyczny sposób rozwiązania tych problemów:

W tym momencie mieliśmy nasz pierwotny milion cyfr, 20 000 kart IBM z 50 cyframi na kartę, z niewielkim, ale zauważalnym nieparzystym wyrównaniem ujawnionym przez analizę statystyczną. Zdecydowano teraz o ponownym losowaniu stołu lub przynajmniej zmianie go przez małą grę w ruletkę, aby usunąć nieparzyste wyrównanie. Dodaliśmy (mod 10) cyfry na każdej karcie, cyfra po cyfrze, do odpowiednich cyfr poprzedniej karty. Otrzymana tabela zawierająca milion cyfr została następnie poddana różnym standardowym testom, testom częstotliwości, testom seryjnym, testom pokera itp. Te miliony cyfr mają czysty stan zdrowia i zostały przyjęte jako nowoczesna tabela losowych cyfr RAND.

Oczywiście istniał dobry powód, by sądzić, że proces dodawania przyniesie jakieś korzyści. Ogólnie rzecz biorąc, mechanizmem leżącym u podstaw jest podejście ograniczające sumy zmiennych losowych modulo interwał jednostkowy w rozkładzie prostokątnym, w taki sam sposób, jak nieograniczone sumy zmiennych losowych zbliżają się do normalności. Metodę tę zastosowali Horton i Smith z Interstate Commerce Commission, aby uzyskać dobre partie pozornie losowych liczb z większych partii złych liczb losowych.

$[0,1]$$u$$[0,1]$$F$$F^{-1}(u)$

Referencje

• Badger, L. (1994). „ Szczęśliwe przybliżenie π Lazzariniego ”. Magazyn Matematyki . Mathematical Association of America. 67 (2): 83–91.

• $(*)$

• Darwin, GH (1877). „ O omylnych pomiarach zmiennych wielkości i traktowaniu obserwacji meteorologicznych. ” Philosophical Magazine , 4 (22), 1–14

• De Forest, EL (1876). Interpolacja i regulacja szeregów . Tuttle, Morehouse and Taylor, New Haven, Conn.

• Galton, F. (1890). „Kości do eksperymentów statystycznych”. Nature , 42 , 13-14

• Stigler, SM (1991). „Symulacja stochastyczna w XIX wieku”. Nauki statystyczne , 6 (1), 89-97.

$(*)$„Każdy, kto rozważa arytmetyczne metody tworzenia losowych cyfr, jest oczywiście w stanie grzechu. Bo, jak już wielokrotnie podkreślano, nie istnieje coś takiego jak liczba losowa - istnieją tylko metody do generowania liczb losowych , a ścisła procedura arytmetyczna oczywiście nie jest taką metodą ”.

Na przykład w moim systemie operacyjnym Linux mogę to sprawdzić

$ date +%s.%N
1479733744.077762986
$ date +%s.%N
1479733980.615056616

> sqrt(-2*log(.077762986))*cos(2*pi*.615056616)
[1] -1.694815

Dodatek: ponieważ obliczanie logarytmów i cosinusów można uznać za niewystarczające manualnie , istnieje wariant Boxa-Müllera, który unika korzystania z tych funkcji transcendentalnych (patrz Ćwiczenie 2.9 w naszej książce Metody statystyczne Monte Carlo ):

Teraz można argumentować przeciwko tej wersji ze względu na zmienne wykładnicze. Istnieje jednak bardzo sprytny sposób symulowania tych zmiennych bez wywołania funkcji transcendentalnych , dzięki von Neumannowi, jak podsumowano w tym algorytmie odtworzonym z niejednorodnych losowych wariantów Luca Devroye'a :

Wprawdzie wymaga to obliczenia 1 / e, ale tylko raz.

$(0,1)^2$$(0,1)$$(0,2\pi)$

Korzystanie z CLT w celu przybliżenia normalności z pewnością nie jest metodą, którą zaleciłbym, ponieważ (1) nadal potrzebujesz innych zmiennych, aby podać średnią, więc równie dobrze możesz użyć mundurów w algorytmie Boxa-Müllera i (2) dokładność rośnie dość powoli z liczbą symulacji. Zwłaszcza jeśli używasz dyskretnej zmiennej losowej, takiej jak wynik kości, nawet z więcej niż sześcioma twarzami . Cytat z Thomas i in. (2007), badanie zalet i wad gaussowskich generatorów losowych:

Centralne twierdzenie graniczne jest oczywiście przykładem metody „przybliżonej” - nawet jeśli zastosowana zostanie idealna arytmetyka, dla skończonego K wynik nie będzie Gaussa.

Oto szybki eksperyment ilustrujący problem: wygenerowałem 100-krotność średniej z 30 wyników kości:

dies=apply(matrix(sample(1:6,30*100,rep=TRUE),ncol=30),1,mean)

następnie znormalizował te średnie do średniej zero - wariancja zmienia się

stdies=(dies-3.5)/sqrt(35/12/30)

i spojrzał na normalne dopasowanie [lub jego brak] tej próbki:

dies76/30122/30$D_i$

dies=matrix(apply(matrix(sample(0:5,15*200,rep=TRUE),nrow=15)/6^(1:15),2,sum),ncol=2) 
norma=sqrt(-2*log(dies[,1]))*c(cos(2*pi*dies[,2]),sin(2*pi*dies[,2]))

dopasowanie jest tak dobre, jak można oczekiwać dla Normalnej próbki o wielkości 200 (wystarczy wykreślić inną dla prawdziwej próbki normalnej norma=rnorm(100)):

jak dodatkowo pokazuje test Kołmogorowa-Smirnowa:

> ks.test(norma,pnorm)

        One-sample Kolmogorov-Smirnov test

data:  norma
D = 0.06439, p-value = 0.3783
alternative hypothesis: two-sided

To nie jest dokładnie przypadkowe, ale powinno być wystarczająco blisko, ponieważ wydaje się, że chcesz trudnego eksperymentu.

Użyj telefonu, aby ustawić chronometr. Po dobrych 10 sekundach przestań (im dłużej czekasz, tym bardziej zbliżasz się do naprawdę „losowego” wyniku, ale 10 sekund jest w porządku). Weź ostatnie cyfry (na przykład 10,67 s da ci 67). Zastosuj tabelę percentyli do rozkładu normalnego. W tym przykładzie musisz po prostu wyszukać 0,67, a znajdziesz numer. W takim przypadku wartość wynosi około 0,45. Nie jest to całkowicie precyzyjne, ale da ci solidne oszacowanie.

Jeśli rzucisz poniżej 50, po prostu zrób 100- [Twój wynik] i skorzystaj z tabeli. Twój wynik będzie taki sam, ze znakiem minus, ze względu na symetrię N (0,1).

Rzućmy obiektywną monetą $n$czasy. Zaczynając od zera, liczymy$+1$ jeśli głowy, $-1$jeśli ogony. Po$n$ Moneta przewraca, dzielimy licznik przez $\sqrt n$. Wykorzystując centralne twierdzenie graniczne , jeśli$n$ jest wystarczająco duży, wówczas powinniśmy mieć „przybliżoną realizację” znormalizowanego gaussowskiego $N (0,1)$.

Czemu? Pozwolić

być iid zmiennymi losowymi Bernoulliego z $\mathbb P (X_k = \pm 1) = \frac 12$. Stąd,

Pozwolić $Y := X_1 + X_2 + \cdots + X_n$. Stąd,

Normalizujący,

otrzymujemy zmienną losową o wariancji jednostkowej

Warto zauważyć, że po wygenerowaniu jednorodnej zmiennej (0,1) można wygenerować dowolną zmienną losową, dla której można obliczyć odwrotną wartość cdf, po prostu podłączając jednolitą zmienną losową do odwrotnego CDF.

Jak więc ręcznie obliczyć mundur (0,1)? Jak wspomniano w @Silverfish, istnieje wiele kości używanych przez tradycyjnych graczy RPG. Jednym z nich jest dziesięciostronna kostka. Zakładając, że to uczciwa kostka, możemy teraz wygenerować dyskretny mundur (0, 9).

Możemy również użyć tego munduru (0,9) do przedstawienia pojedynczej cyfry zmiennej losowej. Więc jeśli użyjemy dwóch kości, otrzymamy jednolitą zmienną losową, która może przyjmować wartości$0.01, 0.02, ..., 0.99, 1.00$. Dzięki trzem kostkom możemy uzyskać jednolity rozkład$0.001, 0.002, ..., 0.999, 1.000$.

Możemy więc zbliżyć się do ciągłego munduru (0,1), przybliżając go dyskretnie rozłożonym dyskretnie jednolitym rozkładem z kilkoma kostkami 10-stronnymi. Można to następnie podłączyć do odwrotnego CDF, aby wygenerować losową zmienną będącą przedmiotem zainteresowania.