Generowanie równomiernie rozmieszczonych liczb losowych za pomocą monety

Masz jedną monetę. Możesz go obrócić tyle razy, ile chcesz.

Chcesz wygenerować losową liczbę taką, że gdzie .a ≤ r < b r , a , b $r$$a \leq r < b$$r,a,b\in \mathbb{Z}^+$

Rozkład liczb powinien być jednolity.

Łatwo jest, jeśli :$b -a = 2^n$

r = a + binary2dec(flip n times write 0 for heads and 1 for tails) 

Co jeśli ?$b-a \neq 2^n$

To, czego szukasz, opiera się na próbkowaniu odrzucenia lub metodzie akceptacji-odrzucenia (pamiętaj, że strona Wiki jest nieco techniczna).

Ta metoda jest przydatna w takich sytuacjach: chcesz wybrać jakiś losowy obiekt ze zbioru ( w twoim przypadku losowa liczba całkowita w zestawie ), ale nie wiesz jak to zrobić, ale może wybrać jakiś losowy obiekt z większego zestawu zawierającego pierwszy zestaw (w twoim przypadku dla niektórych przykład ; odpowiada to rzutom monet).[ a , 2 k + a ] k 2 k + a ≥ b $[a,b]$$[a, 2^k + a]$$k$$2^k + a \ge b$$k$

W takim scenariuszu po prostu wybierasz elementy z większego zestawu, aż losowo wybierzesz element z mniejszego zestawu. Jeśli twój mniejszy zestaw jest wystarczająco duży w porównaniu do większego zestawu (w twoim przypadku zawiera najwyżej dwa razy więcej liczb całkowitych niż , co jest wystarczająco dobre), jest to wydajne.[ a , b $[a, 2^k + a]$$[a,b]$

Alternatywny przykład: załóżmy, że chcesz wybrać losowy punkt w okręgu o promieniu 1. Nie jest tak łatwo znaleźć bezpośrednią metodę tego. Zwracamy się do metody akceptowania-odrzucania: próbkujemy punkty na kwadracie 1x1 obejmującym okrąg i sprawdzamy, czy narysowana liczba leży w okręgu.

(dane techniczne: odpowiedź pasuje do wyboru liczby )$a \le x < b$

Ponieważ możesz rzucić monetą tyle razy, ile zechcesz, możesz uzyskać swoje prawdopodobieństwo jak najbardziej zbliżone do tego, co chcesz ujednolicić, wybierając ułamek (używając binarnego podstawnika: odwracasz monetę dla każdej cyfry po punkcie) i pomnóż przez aby uzyskać liczbę od 0 do [ba-1] (zaokrąglając w dół do liczby całkowitej). Dodaj ten numer do i gotowe.r b - a $r\in [0,1]$$r$$b-a$$a$

Przykład : powiedzmy . 1/3 w systemie dwójkowym to 0,0101010101 .... Zatem, jeśli twój rzut wynosi od 0 do 0,010101 ... twój wybór to . jeśli jest pomiędzy 0,010101 .. a 0,10101010 ... twój wybór będzie równy , w przeciwnym razie będzie .b a + 1 a + $b-a=3$$b$$a+1$$a+2$

Jeśli przerzucić monety razy, a następnie każdą liczbę pomiędzy i będą wybrane z prawdopodobieństwem .a b $t$$a$$b$$\frac{1}{b-a}\pm 2^{-(t+1)}$

Wybierz liczbę z następnej większej potęgi 2 zakresu i odrzuć odpowiedzi większe niż .$b-a$

n = b-a;
N = round_to_next_larger_power_of_2(n)
while (1) {
  x = random(0 included to N excluded);
  if (x < n) break;
}
r = a + x;

$D$$D$$k$$D$$d \in D$$d$$A/2^k$$A$$A/2^k \neq 1/{|D|}$

$n$$D$$n \log_2 |D|$$H$$nH$

Jeśli ba nie jest potęgą 2, być może będziesz musiał przerzucić wiele monet, aby uzyskać wynik. Możesz nawet nigdy nie uzyskać rezultatu, ale jest to mało prawdopodobne.

Metody

Najprostszą metodą jest wygenerowanie liczby w [a, a + 2 ^ n), gdzie 2 ^ n> = ba, dopóki nie wyląduje się w [a, b). Metodą tą wyrzucasz dużo entropii.

Droższa metoda pozwala zachować całą entropię, ale staje się bardzo droga obliczeniowo w miarę wzrostu liczby rzutów monetą / rzutów kostką. Intuicyjnie to tak, jakby traktować monety jako cyfry liczby binarnej po prawej stronie przecinka dziesiętnego, konwertując tę ​​liczbę z bazy 2 na bazę ab po i zwracając cyfry tej liczby, gdy „utknęły”.

Przykład

Poniższy kod konwertuje rzuty rzetelnej n-stronnej kostki na rzuty rzetelnej m-stronnej kostki (w twoim przypadku n = 2, m = ab), przy rosnącym koszcie krańcowym wraz ze wzrostem liczby rzutów. Zwróć uwagę na potrzebę podania liczby wymiernej z dowolną dokładnością. Jedną fajną właściwością jest to, że konwersja z n-stronnego na m-stronny iz powrotem do n-stronnego zwróci oryginalny strumień, choć być może opóźniony o kilka rzutów z powodu zablokowania cyfr.

public static IEnumerable<BigInteger> DigitConversion(this IEnumerable<BigInteger> inputStream, BigInteger modIn, BigInteger modOut) {
    //note: values are implicitly scaled so the first unfixed digit of the output ranges from 0 to 1
    Rational b = 0; //offset of the chosen range
    Rational d = 1; //size of the chosen range
    foreach (var r in inputStream) {
        //narrow the chosen range towards the real value represented by the input
        d /= modIn;
        b += d * r;
        //check for output digits that have become fixed
        while (true) {
            var i1 = (b * modOut).Floor();
            var i2 = ((b + d) * modOut).Floor(); //note: ideally b+d-epsilon, but another iteration makes that correction unnecessary
            if (i1 != i2) break; //digit became fixed?
            //fix the next output digit (rescale the range to make next digit range from 0 to 1)
            d *= modOut;
            b *= modOut;
            b -= i1;
            yield return i1;
        }
    }
}

Wygeneruj binarny dziesiętny. Zamiast jawnego przechowywania, wystarczy śledzić minimalne i maksymalne możliwe wartości. Gdy te wartości znajdą się w tej samej liczbie całkowitej, zwróć tę liczbę całkowitą. Szkic kodu znajduje się poniżej.

(Edytuj) Pełniejsze wyjaśnienie: Powiedz, że chcesz wygenerować losową liczbę całkowitą od 1 do 3 włącznie z prawdopodobieństwem 1/3. Robimy to, generując losową rzeczywistą liczbę binarną dziesiętną x w zakresie (0, 1). Jeśli x <1/3, zwraca 1, w przeciwnym razie jeśli x <2/3 zwraca 2, w przeciwnym razie zwraca 3. Zamiast jawnego generowania cyfr dla x, po prostu śledzimy minimalne i maksymalne możliwe wartości x. Początkowo minimalna wartość x wynosi 0, a maksymalna to 1. Jeśli najpierw odwrócisz głowy, pierwsza cyfra x za przecinkiem dziesiętnym (binarnie) wynosi 1. Minimalna możliwa wartość x (binarnie) następnie wynosi 0,100000 = 1/2, a maksimum to 0.111111111 = 1. Teraz, jeśli twój następny rzut to reszka, x zaczyna się od 0.10. Minimalna możliwa wartość to 0,1000000 = 1/2, a maksymalna to 0,1011111 = 3/4. Minimalna możliwa wartość x wynosi 1/2, więc wiesz, że ” nie ma szansy na zwrócenie 1, ponieważ wymaga to x <1/3. Nadal możesz zwrócić 2, jeśli x kończy się na 1/2 <x <2/3 lub 3, jeśli 2/3 <x <3/4. Załóżmy teraz, że trzeci rzut to ogony. Wtedy x musi zaczynać się od 0.100. Min. = 0,10000000 = 1/2, a maks. = 0,100111111 = 5/8. Odkąd 1/3 <1/2 <5/8 <2/3 wiemy, że x musi mieścić się w przedziale (1/3, 2/3), więc możemy przestać generować cyfry x i po prostu zwrócić 2.

Kod robi to w zasadzie, z wyjątkiem tego, że zamiast generować x między 0 a 1, generuje x między a i b, ale zasada jest taka sama.

def gen(a, b):
  min_possible = a
  max_possible = b

  while True:
    floor_min_possible = floor(min_possible)
    floor_max_possible = floor(max_possible)
    if max_possible.is_integer():
      floor_max_possible -= 1

    if floor_max_possible == floor_min_possible:
      return floor_max_possible

    mid = (min_possible + max_possible)/2
    if coin_flip():
      min_possible = mid
    else:
      max_possible = mid

Uwaga: Przetestowałem ten kod pod kątem metody akceptowania / odrzucania i oba dają jednolite rozkłady. Ten kod wymaga mniejszej liczby rzutów monetą niż akceptacja odrzucenia, z wyjątkiem sytuacji, gdy b - a jest bliskie następnej potędze 2. Np. Jeśli chcesz wygenerować a = 0, b = 62, wtedy lepiej zaakceptować / odrzucić. Udało mi się udowodnić, że ten kod może użyć średnio nie więcej niż 2 rzutów monetą niż zaakceptować / odrzucić. Z mojego czytania wynika, że ​​Knuth i Yao (1976) podali metodę rozwiązania tego problemu i udowodnili, że ich metoda jest optymalna pod względem oczekiwanej liczby rzutów monetą. Udowodnili ponadto, że oczekiwana liczba rzutów musi być większa niż entropia rozkładu Shannona. Nie mogłem jednak znaleźć kopii tekstu artykułu i byłbym ciekawy, jaka jest ich metoda. (Aktualizacja: właśnie znalazłem wystawę Knuth Yao 1976 tutaj:http://www.nrbook.com/devroye/Devroye_files/chapter_fifteen_1.pdf, ale jeszcze tego nie przeczytałem). Ktoś wspomniał także o Han Hoshi w tym wątku, który wydaje się bardziej ogólny, i rozwiązuje go za pomocą stronniczej monety. Zobacz także http://paper.ijcsns.org/07_book/200909/20090930.pdf autor: Pae (2009), aby uzyskać dobre omówienie literatury.

Prosta odpowiedź?

$b-a$$r$

$r$$b - a = 2^n + 1$$n$

Jest to proponowane rozwiązanie w przypadku, gdy b - a nie jest równe 2 ^ k. Ma działać w określonej liczbie kroków (nie trzeba wyrzucać kandydatów, którzy znajdują się poza oczekiwanym zakresem).

Nie jestem jednak pewien, czy to poprawne. Prosimy o krytykę i pomoc w opisaniu dokładnej niejednorodności w tym generatorze liczb losowych (jeśli istnieje) oraz w jaki sposób ją zmierzyć / skwantyfikować.

Najpierw przekonwertuj na równoważny problem generowania równomiernie rozmieszczonych liczb losowych w zakresie [0, z-1], gdzie z = b - a.

Niech m = 2 ^ k będzie najmniejszą potęgą 2> = z.

Zgodnie z powyższym rozwiązaniem mamy już równomiernie rozłożony generator liczb losowych R (m) w zakresie [0, m-1] (można to zrobić, podrzucając k monet, po jednej na każdy bit).

    Keep a random seed s and initialize with s = R(m).   

    function random [0, z-1] :
        x = R(m) + s 
        while x >= z:
            x -= z
        s = x
        return x

Pętla while działa najwyżej 3 razy, podając następną liczbę losową w ustalonej liczbie kroków (najlepszy przypadek = najgorszy przypadek).

Zobacz program testowy dla liczb [0,2] tutaj: http://pastebin.com/zuDD2V6H

Teoretycznie optymalny algorytm

Oto poprawka innej odpowiedzi, którą opublikowałem. Druga odpowiedź ma tę zaletę, że łatwiej jest rozszerzyć ją na bardziej ogólny przypadek generowania jednego dyskretnego rozkładu z drugiego. W rzeczywistości druga odpowiedź to szczególny przypadek algorytmu spowodowany przez Hana i Hoshi.

Algorytm, który opiszę tutaj, oparty jest na Knuth i Yao (1976). W swoim artykule udowodnili również, że ten algorytm osiąga minimalną możliwą oczekiwaną liczbę rzutów monetą.

Aby to zilustrować, rozważ metodę próbkowania odrzucenia opisaną przez inne odpowiedzi. Jako przykład załóżmy, że chcesz równomiernie wygenerować jedną z 5 liczb [0, 4]. Następna potęga 2 to 8, więc rzucisz monetą 3 razy i wygenerujesz losową liczbę do 8. Jeśli liczba wynosi od 0 do 4, to ją zwrócisz. W przeciwnym razie wyrzucasz go i generujesz kolejną liczbę do 8 i próbujesz ponownie, aż ci się uda. Ale kiedy wyrzucasz numer, po prostu zmarnowałeś trochę entropii. Zamiast tego możesz uzależnić się od liczby, którą wyrzuciłeś, aby zmniejszyć liczbę przyszłych rzutów monetami, których będziesz oczekiwać. Konkretnie, po wygenerowaniu liczby [0, 7], jeśli jest to [0, 4], zwróć. W przeciwnym razie jest to 5, 6 lub 7 i w każdym przypadku robisz coś innego. Jeśli jest to 5, ponownie rzuć monetą i zwróć 0 lub 1 na podstawie rzutu. Jeśli to 6, rzuć monetą i zwróć 2 lub 3. Jeśli jest to 7, rzuć monetą; jeśli to głowy, zwróć 4, jeśli ogony zaczną od początku.

Pozostała entropia z naszej pierwszej nieudanej próby dała nam 3 przypadki (5, 6 lub 7). Jeśli po prostu to wyrzucimy, wyrzucamy monety log2 (3). Zamiast tego zachowujemy go i łączymy z wynikiem kolejnego przerzucenia, aby wygenerować 6 możliwych przypadków (5H, 5T, 6H, 6T, 7H, 7T), które od razu spróbujmy ponownie wygenerować ostateczną odpowiedź z prawdopodobieństwem sukcesu 5/6 .

Oto kod:

# returns an int from [0, b)
def __gen(b):
  rand_num = 0
  num_choices = 1

  while True:
    num_choices *= 2
    rand_num *= 2
    if coin.flip():
      rand_num += 1

    if num_choices >= b:
      if rand_num < b:
        return rand_num
      num_choices -= b
      rand_num -= b

# returns an int from [a, b)
def gen(a, b):
  return a + __gen(b - a)