Jeśli ktoś grał przykładowo w Diablo (lub inne gry), to być może spotkał się z tabelkami opisującymi z jakim prawdopodobieństwem wypadnie dany przedmiot. Ja sam bawiłem się kiedyś pisaniem gierek (nawet napisałem taką jedną gierkę RPG (totalnie biedna) - ale to inna historia) i spotkałem się z tym właśnie problemem - w jaki sposób skutecznie losować jakąś liczbę(*) z zadanym prawdopodobieństwem?
(*) - Wszystko w informatyce można sprowadzić do liczb. W tym przypadku - indeks w tablicy z przedmiotami.
Podejście naiwne jest takie, żeby wypisać w tabeli względne częstości każdej rzeczy, następnie wylosować jakąś liczbę (z zakresu od 0 do suma_częstości-1) i odejmować każde z tych prawdopodobieństw aż liczba zrobi się mniejsza od 0. Wtedy dana rzecz zostaje wylosowana.
Można tę metodę zaimplementować niejako w drugą stronę: tabelę z częstościami przekształcić w taki sposób, aby dla każdej rzeczy zsumować jej częstość z częstościami wszystkich poprzednich rzeczy. Następnie wylosować jakąś liczbę i porównywać po kolei z sumowanymi częstościami - jeśli liczba zrobi się większa lub równa danej sumie, losujemy poprzednią rzecz.
Podejście naiwne
[spis treści]
int n; //dana ilość rzeczy
int prob[n]; //dana tablica ze względnymi częstościami
int losuj[n]; //tworzona tablica z sumowanymi częstościami
//inicjacja tablicy losującej
void inicjuj()
{
int i;
for(i=1; i<n; ++i) losuj[i]=prob[i-1]+prob[i];
}
//funkcja losująca
int losuj_pozycje()
{
int i, p;
p=random(losuj[n-1]);
for(i=0; losuj[i]<p; ++i);
return i;
}
Podejście naiwne jest bardzo łatwe do zrealizowania, ale jest też wolne - ma złożoność liniową. Dla dużej ilości pozycji i dla dużej ilości losowań spowolnienie wykonywania programu może być zauważalne.
Niedawno jednak natrafiłem na algorytm, który po zbudowaniu (dosyć szybkim) tabeli losującej pozwala na losowanie danej pozycji ze złożonością O(1) - innymi słowy jest tak szybkie, jak zwykłe losowanie. Jako ciekawostkę powiem, że na ten algorytm natrafiłem w artykule dotyczącym ... genetyki. Jak widać informatyka nie znajduje zastosowania tylko w grach .
Algorytm ten nazywa się „Metoda Aliasów Vosego” i należy do rodziny metod aliasów pozwalających szybko losować wg. dyskretnego rozkładu prawdopodobieństwa (inaczej mówiąc - wg. zadanego dla każdego elementu prawdopodobieństwa).
Opis metody
[spis treści]Metoda jest (na ile mogę to stwierdzić) modyfikacją metody Monte Carlo. Można to sobie wyobrazić jako rzucanie strzałkami do tarczy, na której narysowane są słupki o wysokości odpowiadającej poszczególnym prawdopodobieństwom. W każdy słupek trafimy zgodnie z odpowiednim prawdopodobieństwem, jednak czasem trafimy na puste pole - musimy powtarzać wtedy rzut. Aby tego uniknąć, równamy wszystkie słupki do średniej w taki sposób, że od wystających słupków odcinamy odpowiednio duże kawałki, którymi wypełniamy brakujące przestrzenie nad małymi słupkami. Nad każdym słupkiem może być co najwyżej jeden alias.
|
|
| Ilustracja metody Monte Carlo. Po trafieniu w puste miejsce (szare pole) trzeba powtórzyć losowanie, co jest czasochłonne. | Ilustracja metody aliasów. Tablica została zmieniona tak, że nadmiarowe słupki wypełniają wolne miejsca nad słupkami niedostającymi do średniej wysokości. |
Budowa tablicy
[spis treści]Algorytm wzięty z artykułu Keitha Schwartza Darts, Dice, and Coins: Sampling from a Discrete Distribution. Finalną wersję oraz implementację algorytmu w Javie można znaleźć na końcu artykułu.
- Stwórz tablice Alias i Prob, każda rozmiaru n.
- Stwórz dwie listy Small i Large(*).
- Wymnóż każde prawdopodobieństwo przez n(**).
-
Dla każdego przeskalowanego prawdopodobieństwa p[i]:
- Jeśli p[i]<1, dodaj element p[i] do listy Small.
- W przeciwnym wypadku (p[i]≥1) do listy Large.
-
Dopóki listy Small oraz Large zawierają jakieś elementy:
- Wyciągnij z listy Small pierwszy element (do zmiennej l).
- Wyciągnij z listy Large pierwszy element (do zmiennej g).
- Ustaw Prob[l]=p[l].
- Ustaw Alias[l]=g.
- Ustaw p[g]=(p[g]+p[l])−1.
- Jeśli p[g]<1, dodaj g do listy Small.
- W przeciwnym wypadku (p[g]≥1), dodaj g do listy Large.
-
Dopóki Large jest niepusta:
- Wyciągnij pierwszy element z listy Large (do zmiennej g).
- Ustaw Prob[g]=1.
-
Dopóki Small jest niepusta(***):
- Wyciągnij pierwszy element z listy Small (do zmiennej l).
- Ustaw Prob[l]=1.
(*) - Listy można w bardzo łatwy sposób zaimplementować tablicami:
int listaSize=0, Lista[n]; //dodawanie Lista[listaSize++]=x; //usuwanie x=Lista[--listaSize];
Przy odrobinie inwencji dwie listy Small i Large można by zaimplementować na jednej tablicy, ale uważam, że jest to zbytek komplikowania sobie życia. Już lepiej użyć list / kolejek.
(**) - Wymnażanie prawdopodobieństw ma na celu ich znormalizowanie - tak, aby średnia wartość znormalizowanych prawdopodobieństw była równa 1. Suma pierwotnych (wejściowych) prawdopodobieństw musi się równać 1. Algorytm można jednak przystosować dla dowolnych liczb - także całkowitych (dla względnych częstości).
(***) - Wg. Keitha lista Small może zostać niepusta tylko w skutek numerycznej niestabilności (chodzi o dokładność liczb zmiennoprzecinkowych).
Punkty 1-4 to inicjacja danych. Punkt 5 to główna pętla budująca tablicę aliasów. Punkty 6 i 7 mają za zadanie ustawienie tych elementów, które z jakiś powodów nie brały udziału w budowaniu tablicy. Normalnie powinno to być możliwe tylko dla elementów dla których prawdopodobieństwo jest równych średniej (p[g]=1).
Losowanie liczby przy pomocy tej metody polega na tym, aby najpierw wylosować pozycję. Gdy już mamy daną pozycję, sprawdzamy za pomocą rzutu fałszywą monetą, czy aby na pewno chodzi nam o tę pozycję, czy może o jej alias.
Losowanie
[spis treści]- Losuj liczbę całkowitą i od 0 do n-1.
- Losuj liczbę rzeczywistą r od 0 do 1 i porównaj z prawdopodobieństwem Prob[i].
- Jeśli r < Prob[i], zwróć i.
- W przeciwnym wypadku zwróć Alias[i].
Implementacja
[spis treści]Implementacja tego algorytmu jest naprawdę prosta. Pozwolę sobie przedstawić dwie implementacje (używające tablic jako list). Pierwsza - dla prawdopodobieństw zmiennoprzecinkowych (tzw. - dowolnych), druga - dla prawdopodobieństw całkowitych (częstości względnych).
Wersja zmiennoprzecinkowa
[spis treści]
int Alias[];
double Prob[];
int n=0;
/*
Funkcja obliczająca tablice aliasów.
Oczekuje jako parametru tablicy (o wielkości n)
z prawdopodobieństwami.
Mogą być dowolne (nie sumować się do 1) -
do kodu dodałem kawałek, który je normalizuje.
*/
void ObliczTablice(int size, double dListaLiczb[])
{
int nSmall, Small[]; //lista Small
int nLarge, Large[]; //lista Large
double p[]; //znormalizowane prawdopodobieństwa
double dSuma;
int i, l, g;
//1,2 - inicjalizacja tablic/list
n=size;
p=new double[n];
Small=new int[n];
nSmall=0;
Large=new int[n];
nLarge=0;
//3 - zsumowanie prawdopodobieństw - w celu ich normalizacji
for(i=0, dSuma=0; i<n; ++i) {
p[i]=dListaLiczb[i];
dSuma+=p[i];
}
//4 - normalizacja i wypełnienie list Small i Large
for(i=0; i<n; ++i) {
p[i]*=(double)n; //z algorytmu
p[i]/=dSuma; //normalizacja do sumy prawdopodobieństw=1
if(p[i]<1) Small[nSmall++]=i;
else Large[nLarge++]=i;
}
//5 - budowanie listy Alias i Prob
while(nSmall>0 && nLarge>0) {
l=Small[--nSmall];
g=Large[--nLarge];
Prob[l]=p[l];
Alias[l]=g;
p[g]=(p[g]+p[l])-1;
if(p[g]<1) Small[nSmall++]=g;
else Large[nLarge++]=g;
}
//6,7 - końcowe sprzątanie
while(nLarge>0) {
g=Large[--nLarge];
Prob[g]=1;
}
while(nSmall>0) {
l=Small[--nSmall];
Prob[l]=1;
}
}
/*
Funkcja losująca.
rand()/RAND_MAX można zastąpić dowolną funkcją zwracającą
liczbę zmiennoprzecinkową z zakresu [0.0, 1.0)
*/
int LosujPozycje(void)
{
int i;
double p;
if(n==0) return 0; //zabezpieczenie, gdy tablica
//losująca niezainicjowana
i=random(n); //zwraca liczbę całkowitą z zakresu 0÷(n-1)
p=(double)rand()/(double)RAND_MAX;
if(p<Prob[i]) return i;
else return Alias[i];
}
Wersja ta operuje na liczbach zmiennoprzecinkowych, co wiąże się z ryzykiem występowania niedokładności (niestabilności algorytmu). Można ją w łatwy sposób przerobić na wersję całkowitoliczbową. Zamiast porównywać z 1 trzeba będzie porównywać częstości z sumą (konkretnie: ze średnią pomnożoną przez n - czyli sumą właśnie).
Wersja całkowitoliczbowa
[spis treści]
int Alias[];
int Prob[];
int n=0, s; //s - średnia pomnożona przez n, czyli suma
/*
Funkcja obliczająca tablice aliasów.
Oczekuje jako parametru tablicy
z częstościami względnymi podanymi jako liczby całkowite.
*/
void InicjujTablice(int size,int iListaLiczb[])
{
int nSmall, Small[]; //lista Small
int nLarge, Large[]; //lista Large
int p[]; //znormalizowane prawdopodobieństwa
int i, l, g;
//1,2 - inicjalizacja tablic/list
n=size;
Alias=new int[n];
Prob=new int[n];
p=new int[n];
Small=new int[n];
nSmall=0;
Large=new int[n];
nLarge=0;
//3 - zsumowanie prawdopodobieństw - w celu ich normalizacji
for(i=0, s=0; i<n; ++i) {
p[i]=iListaLiczb[i];
s+=p[i];
}
//4 - normalizacja i wypełnienie list Small i Large
for(i=0; i<n; ++i) {
p[i]*=n; //z algorytmu
//p[i]/=s; //nie zachodzi, patrz przypis(*)
if(p[i]<s) Small[nSmall++]=i;
else if(p[i]>s) Large[nLarge++]=i;
else Prob[i]=s; //patrz przypis(**)
}
//5 - budowanie listy Alias i Prob
while(nSmall>0 && nLarge>0) {
l=Small[--nSmall];
g=Large[--nLarge];
Prob[l]=p[l];
Alias[l]=g;
p[g]=(p[g]+p[l])-s;
if(p[g]<s) Small[nSmall++]=g;
else if(p[g]>s) Large[nLarge++]=g;
else Prob[g]=s; //jak wyżej(**)
}
//6,7 - końcowe sprzątanie - nie wiem, czy potrzebne
while(nLarge>0) {
g=Large[--nLarge];
Prob[g]=s;
}
while(nSmall>0) {
l=Small[--nSmall];
Prob[l]=s;
}
}
/*
Funkcja losująca.
*/
int LosujPozycje(void)
{
int i, p;
if(n==0) return 0; //zabezpieczenie, gdy tablica
//losująca niezainicjowana
i=random(n); //zwraca liczbę całkowitą z zakresu 0÷(n-1)
p=random(s);
if(p<Prob[i]) return i;
return Alias[i];
}
(*) - Różnica między tą wersją a zmiennoprzecinkową jest taka, że nie dzieli się prawdopodobieństw przez sumę - dzięki temu możemy operować na całkowitych liczbach. Można by pokusić się o podzielenie wszystkich liczb przez NWD żeby nie ryzykować użycia zbyt dużych liczb całkowitych. Ale obliczenie NWD dla wielu liczb może być zbyt czasochłonne i może de facto nie zmienić zbyt wiele, lub nawet niczego (NWD może wyjść 1).
(**) - Druga różnica jest taka, że używając liczb całkowitych możemy mieć pewność, gdy jakiś słupek wypełnia daną pozycję całkowicie - możemy więc od razu ustawić dla niego wartość tablicy Prob. Wartość tablicy Alias może być dowolna - i tak nie będzie w takim przypadku używana.
Istnieje możliwość pozbycia się list Small i Large za pomocą trzech indeksów. Oszczędza to pamięć, chociaż szczerze mówiąc niewiele. Wersję taką znalazłem na Panda's Thumb. Zawiera ona błąd, ale nie mogę znaleźć rozwiązania tego problemu. Wysłałem co prawda maila do autora artykułu i chociaż poprawił mniejsze pomyłki i zapewnił mnie, że zbada ów błąd, to jak do tej pory nic się w samym algorytmie nie zmieniło.
Pozwolę sobie przepisać tę wersję ze strony angielskiej (z drobnymi zmianami).
Wersja bez listy Small i Large (zawiera błąd)
[spis treści]
int a[]; //aliasy
long h[]; //”wysokości słupków” - 64-bitowy int, patrz przypis(*)
int n=0;
/*
* Zainicjuj tablicę aliasów z zadanej dystrybucji dyskretnej, pp
*/
void create_alias_table(int size, const double *pp) {
double f=0.0, p[];
int i,g,m,mm;
n=size;
a=new int[n];
p=new int[n];
// znormalizuj pp
for(i=0;i<n;++i)
f += pp[i];
f = (double)n/f;
for(i=0;i<n;++i)
p[i] = pp[i]*f;
// znajdź pozycje startowe
for(g=0; g<n && p[g] < 1.0; ++g)
/*noop*/;
for(m=0; m<n && p[m] >= 1.0; ++m)
/*noop*/;
mm = m+1;
//buduj tablicę aliasów, dopóki nie skończą się małe i duże słupki
while(g < n && m < n) {
// skonwertuj double na 64-bit int, control for precision
h[m] = (long)(ceil(p[m]*9007199254740992.0)) << 11);
a[m] = g;
p[g] = (p[g]+p[m])-1.0;
if(p[g] >= 1.0 || mm < g) {
for(m=mm;m<n && p[m] >= 1.0; ++m)
/*noop*/;
mm = m+1; //patrz przypis(**)
} else {
m = g;
}
//tutaj alternatywna wersja - patrz przypis(**)
for(; g<n && p[g] < 1.0; ++g)
/*noop*/;
}
// jeśli jakieś słupki zostały, nie mają aliasów
for(; g<n; ++g) {
if(p[g] < 1.0)
continue;
h[g] = MAX_LONG;
a[g] = g;
}
if(m < n) {
h[m] = MAX_LONG;
a[m] = m;
for(m=mm; m<n; ++m) {
if(p[m] > 1.0)
continue;
h[g] = MAX_LONG;
a[m] = m;
}
}
}
/*
Funkcja losująca.
Autor używa funkcji losującej skonstruowanej specjalnie
do konkretnego zadania. Zmieniłem ją, nie wiem jednak,
czy ta zmiana nie wpłynie na jakość jej losowości.
*/
int get_position() {
long u = rand_uint64();
int x = u%64; //dokonana przeze mnie zmiana,
//można użyć po prostu kolejnego losowania
return (u < h[x]) ? x : a[x];
}
(*) - Autor używa dużych, 64-bitowych liczb całkowitych dla oddania wysokości słupków, tj. prawdopodobieństwa. Robi to ze względu na precyzję, ale myślę, że wyrzucenie tablicy h[] i zastąpienie jej tablicą p[] w zupełności spełni swoje zadanie do zwykłych zastosowań (maksymalną wartością w p[] powinna być 1.0).
(**) - Gdy próbowałem zrozumieć działanie tego algorytmu, zauważyłem, że dla pewnych danych działa błędnie i w związku z tym produkuje błędne tablice aliasów. Miało to związek z wykonywaniem działania mm=m+1 tylko dla jednego przypadku. Wyrzucenie tego poza if-else poprawiało działanie dla tamtego zestawu danych, ale za to powodowało błędne obliczenia dla innych zestawów. W tej chwili nie mam pojęcia, jak poprawić ten kod.
Końcowe uwagi
[spis treści]W zależności od tego, który element będziemy wyciągać z listy Small i Large (ostatni / pierwszy / dowolny) końcowe tablice Aliasów i Prawdopodobieństw mogą się różnić, nie powinno to mieć jednak żadnego znaczenia - tj. efekty funkcji losującej będą takie same.
Starałem się, aby umieszczony kod był w pełni wykonywalny, chociaż miejscami zmieniałem kod już w tekście - przez co czasem mógł się zawieruszyć tu i ówdzie błąd (literówka, czy coś). Użyłem języka C z ukłonem w stronę Javy (której się aktualnie uczę). Chodzi głównie o traktowanie zmiennych tablicowych, konkretne właściwości funkcji random, typu long i chyba jeszcze paru rzeczy.
Jednak kopiując podany kod, niezależnie od użytego języka, trzeba będzie i tak dokonać odpowiednich poprawek, żeby współgrał on z resztą programu. Liczę na to, że osoba znająca podstawy programowania nie będzie miała problemu ze zrozumieniem niuansów zawartych w kodzie oraz z odpowiednią jego modyfikacją.
Linki
Brak komentarzy:
Prześlij komentarz