Generator liczb losowych: rozkład Poissona

Niedawno pisałem o rzeczywistym przykładzie rozkładu Poissona. Dzisiaj chciałbym napisać o tym, w jaki sposób generować programowo liczby losowe o takim rozkładzie.

Na matematyce nie znam się zbyt dobrze, a na kombinatoryce i prawdopodobieństwie prawie w ogóle(*) - z tego powodu nie byłem w stanie wymyślić odpowiedniego algorytmu samemu. Postanowiłem więc poszukać odpowiedniego w internecie. Ale szczerze mówiąc nawet gdybym był w stanie coś samemu wymyślić i tak szukałbym gotowego rozwiązania, ponieważ zdolność do szukania gotowych rozwiązań to jedna z cnót programisty (mówię poważnie).

(*) - Ostatnio jednak poczytałem co nieco o tym rozkładzie i to, co zrozumiałem chciałbym przedstawić kiedyś indziej.

Ten artykuł nie będzie więc zbyt twórczy, ale to nawet lepiej . Przedstawię tu tylko gotowe implementacje w Javie kilku funkcji, które znalazłem na stronie Wikipedii.

Algorytm Knutha

Algorytm ten bazuje na tzw. pseudolosowym samplowaniu. Metody jeszcze nie zrozumiałem , doszedłem tylko do tego, że służą one do generowania pseudolosowych liczb o zadanej dystrybucji oraz że jednym z najprostszych przykładów tego rodzaju algorytmów jest metoda Monte-Carlo, którą akurat rozumiem .

    public int gen_poiss_knuth(int lbd) {
        double L, p, u;
        int k;
        
        L=Math.exp(-lbd);
        k=0;
        p=1;
        do {
            k=k+1;
            u=Math.random();
            p=p*u;
        } while(p>L);
        
        return (k-1);
    }

Nie zmieniałem działań na typowo C-owskie, takie jak ++ czy *=

Jak widać, algorytm jest w miarę prosty i krótki. Ma złożoność liniową względem k, które średnio wynosi lbd.

Algorytm Junhao, oparty na algorytmie Knutha

Powyższy algorytm jest liniowy, co nie jest zbyt optymalne. Istnieją modyfikacje, takie jak poniższa, które są szybsze, ale mogą być niestabilne numerycznie dla dużych wartości lambdy.

    public int gen_poiss_Junhao(int lbd) {
        double lbdLeft, p, u;
        int k;
        double STEP=500;
        
        lbdLeft=lbd;
        k=0;
        p=1;
        do {
            k=k+1;
            u=Math.random();
            p=p*u;
            if(p<Math.E && lbdLeft>0) {
                if(lbdLeft>STEP) {
                    p=p*Math.exp(STEP);
                    lbdLeft=lbdLeft-STEP;
                } else {
                    p=p*Math.exp(lbdLeft);
                    lbdLeft=-1;
                }
            }
        } while(p>1);
        
        return (k-1);
    }

Wartość STEP zależy od długości słowa reprezentującego liczbę zmiennoprzecinkową. Dla typu double 500 powinno być dobrym wyborem (patrz więcej na Wikipedii).

Algorytm odwrotnej transformaty

Algorytm ten jest prosty i efektywny dla niewielkich wartości lambda. Gdy jednak lambda rośnie, jego efektywność staje się proporcjonalna do jej wielkości i dodatkowo zaczynają się kumulować błędy zaokrągleń.

    public int gen_poiss_ITS(int lbd) {
        double p, s, u;
        int x;
        
        x=0;
        p=Math.exp(-lbd);
        s=p;
        u=Math.random();
        while(u>s) {
            x=x+1;
            p=p*lbd/x;
            s=s+p;            
        }
        
        return x;
    }

I to w zasadzie wszystkie algorytmy, które chciałbym przedstawić. Jest jeszcze jeden, który znalazłem, bazujący na naiwnym zrozumieniu problemu, który chciałbym pokazać dla porównania oraz dla ilustracji, dlaczego lepiej szukać gotowych rozwiązań (algorytm poniższy prawdopodobnie byłbym w stanie wymyślić sam gdybym był zrozumiał rozkład Poissona).

Algorytm naiwny

Algorytm ten znalazłem na takiej oto stronie (nie wiem, co to za projekt, ale póki co nie obchodzi mnie to ).

    public int gen_poiss(int mi, int NT) {
        Random rnd=new Random();
        int cnt, r, i;

        cnt=0;
        for(i=0; i<INT; ++i) {
            r=rnd.nextInt(NT);
            if(r<mi) cnt++; //operator '<' albo '<='
                            //patrz porównanie niżej
        }

        return cnt;
    }

Jest to zmodyfikowana wersja - funkcja ta losuje jedną zmienną losową, podczas gdy na stronie przedstawiona jest funkcja wykonująca symulację wielu losowań. Jest ona naiwna, bazująca na bezpośrednim „działaniu” rozkładu Poissona. Z tego powodu nie daje ona zbyt dobrych rezultatów.

Czerwony i zielony słupek – metoda naiwna z operatorem '<=' lub '<', patrz miejsce w kodzie powyżej.
Obrazek wygenerowany w OpenOffice Calc i obrobiony w GIMPie, przeze mnie.

Jak widać na powyższym obrazku, metody znalezione na Wikipedii (niebieskie słupki) pasują praktycznie idealnie do rozkładu Poissona (linia ciągła), podczas gdy ostatnia, naiwna metoda znacząco od niego odbiega.

I to wszystko. Życzę miłego programowania .

PS. Uwaga dla każdego, kto chce umieszczać kod na stronie – nie zapomnijcie zamieniać znaków '<' i na wszelki wypadek również '>' na '&lt;' oraz '&gt;' .