Generowanie liczb losowych

[Crizz]

Mając daną procedurę \(\displaystyle{ RANDOM(0,1)}\), która zwraca 0 lub 1, każde z nich z prawdopodobieństwem \(\displaystyle{ \frac{1}{2}}\), wykorzystać ją do stworzenia procedury \(\displaystyle{ RANDOM(a,b)}\), która zwraca losową liczbę całkowitą z przedziału \(\displaystyle{ <a,b>}\).

Rozkład prawdopodobieństwa ma być jednostajny, tzn. każda z liczb ma być zwracana z prawdopodobieństwem \(\displaystyle{ \frac{1}{b-a+1}}\).

Ma ktoś pomysł?

[Elminster]

Proponuję następującą funkcję:

Najpierw robimy pętlę, która \(\displaystyle{ b-a}\) raza wykonuje funkcję random i wyniki zsumowujemy do jednej zmiennej. Następnie do tej zmiennej dodajemy jeszcze wartość a, i w ten sposób otrzymamy liczbę zawierającą się w przedziale od a do b, a prawdopodobieństwo każdego wyniku będzie takie same (jednostajne)

[Crizz]

prawdopodobieństwo każdego wyniku będzie takie same (jednostajne)

Otóż w tym właśnie problem, że nie. Taki algorytm wykonuje \(\displaystyle{ b-a}\) niezależnych prób losowych i zlicza, ile było wśród nich sukcesów (wylosowań jedynki), przy czym w każdej próbie prawdopodobieństwo sukcesu jest takie same. A to jest typowy rozkład dwumianowy (czyli np. dużo większe prawdopodobieństwo wylosowania przypada liczbie \(\displaystyle{ \frac{a+b}{2}}\) niż liczbie \(\displaystyle{ b}\)).

[drunkard]

Można spróbować tak: dzielisz przedział na dwie części i w zależności od wyniku RANDOM01 zawężasz przedział o połowę, tzn. mniej więcej tak:

niech pocz = a
niech koniec = b
dobre parę razy wykonaj:
if (RANDOM01 równa się 0) to
koniec = (pocz+koniec)/2
w przeciwnym razie
pocz = (pocz+koniec)/2
od
wynikiem jest (pocz+koniec)/2

Rzecz jasna będzie to jakoś tam przybliżony rozkład jednostajny.

[Crizz]

drunkard, pomysł, który podał Elminster przyszedł mi do głowy jako pierwszy, natomiast twój jako drugi .

Niestety, zadanie wyraźnie mówi o dokładnym rozkładzie jednostajnym, a ten algorytm można zbudować tak, żeby działał dokładnie tylko w przypadku, gdy długość przedziału jest potęgą dwójki (przynajmniej tak mi się wydaje).

[drunkard]

Wydaje mi się, że dokładnie to nigdy nie będziesz miał... Im więcej iteracji, tym większa dokładność.. ograniczona zresztą typem danych, jaki wybierzesz na reprezentację liczby rzeczywistej...

[smiechowiec]

Proponowane testy warto sprawdzić dla liczb 1, 2, 3
pokazując, że każda liczba jest wybierana z równym prawdopodobieństwem.

Można zaproponować metodę eliminacji
1. Tworzymy zbiór liczb akceptowanych od a b ( np w tablicy).
2. Dla każdej liczby z aktualnego zbioru wywołujemy test który zwraca liczbę 0 lub 1.
3. Jeśli liczba wyrazów, dla których wypadło 1 jest równa zero wracamy do punku 2
4. Eliminujemy ze zbioru liczby dla których wynik testu wynosi 0.
5. Jeśli pozostała liczba wyrazów jest większa od 1 wracamy do punku 2.
6. Pozostała liczba jest wybraną losowo liczbą bez preferencji jakiegokolwiek zakresu.

[Fibik]

Normalnie losujemy n bitów, np. 31 sztuk z liczby int32, i mamy już losową liczbę r z przedziału <0, range>; range = 2^(n+1)-1; (można też tak zrobić typ real/float <0,1>).

Teraz tylko skalujemy to do przedziału <a,b>:
x = a + (b-a) * r / range;

Albo tak: w = b-a + 1; i tu liczymy bity + 1 (albo mamy stały limit np. n = 32);
Losujemy te n bitów, i to jest losowa liczba r;
x = a + r mod w;

-- 30 grudnia 2009, 16:59 --

tyle że to nie będzie rozkład jednostajny - chociaż prawdop. wylosowania każdej liczby będzie jednakowe: p = 1/N, N - ilość liczb w przedziale;
losujemy ciąg n-1 bitów, więc prawdopodobieństwo otrzymania dowolnej liczby spośród N = 2^n (od 0 do 2^n - 1) wynosi: 1/2 * 1/2 ... = 1/2^n = 1/N.

Natomiast prawdopodobieństwo wylosowania liczby z danego przedziału: p(a < x < b) nie będzie równe (b-a)/N, a tak jest w rozkładzie jednostajnym.

p(x < N/2); górna połowa bitów ma być wyzerowana, a dolna nie ma znaczenia: p = 1/2^(n/2), co będzie raczej różne od N/2/N = 1/2.

Trzeba wstawiać bity tą 'metodą binarną', wtedy wyjdzie jednostajnie.

-- 31 grudnia 2009, 18:43 --

dobra, jeszcze raz:
p(x < N/2) - tylko najstarszy bit ma być 0, zatem p = N/2/N = 1/2, czyli jednak pasuje...
górna połowa bitów zero: x < 2^(n/2) = N^0.5, wtedy p = 1/2^(n/2) = 1/N^0.5 = N^0.5/N = ...