Liczby losowe – metoda odwracania dystrybuanty

06.01.2021

Gęstość prawdopodobieństwa i dystrybuanta

Rozważmy funkcję $f(x)$ określoną na $\mathbb{R}$, nieujemną i całkowalną. Powiemy, że $f(x)$ jest gęstością rozkładu prawdopodobieństwa $P$, jeżeli dla dowolnego zbioru $A\subset \mathbb{R}$ jest

$P(A)=\underset{A}{\int }f(x)dx$

gdzie $P(A)$ jest prawdopodobieństwem przypisanym zbiorowi $A$.

Z powyższej definicji wynika natychmiast, że funkcja $f(x)$ musi być unormowana:

${\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}f(x)dx=1$

Dystrybuantą gęstości prawdopodobieństwa $f(x)$ nazywamy funkcję:

$F(x)={\int }_{-\mathrm{\infty }}^{x}f(t)dt$

Prawdopodobieństwo, że zmienna losowa przyjmie wartość należącą do przedziału $[a,b]$ wynosi:

$P(a⩽x⩽b)={\int }_a^{b}f(t)dt=F(b)-F(a)$

Rozkład równomierny $U(a,b)$

Rozkładem równomiernym (prostokątnym) nazywamy rozkład prawdopodobieństwa, dla którego gęstość prawdopodobieństwa $f(x)$ jest równa:

$f(x)=\{\begin{array}{lll}0& \text{dla}& x<a\\ \frac{1}{b-a}& \text{dla}& x\in [a,b]\\ 0& \text{dla}& x>b\end{array}$

Tak zdefiniowany rozkład równomierny będziemy oznaczali $U(a,b)$. Dystrybuanta tego rozkładu jest określona wzorem:

$F(x)=\{\begin{array}{lll}0& \text{dla}& x<a\\ \frac{x-a}{b-a}& \text{dla}& x\in [a,b]\\ 1& \text{dla}& x>b\end{array}$

Zbiór liczb należących do rozkładu równomiernego $U(0,1)$ możemy łatwo uzyskać. Przykładowo w arkuszu kalkulacyjnym LibreOffice dostępna jest funkcja $\text{LOS}()$, która zwraca przypadkową liczbę z przedziału $[0,1)$. Pisząc makro mamy dostępną analogiczną funkcję $\text{Rnd}()$. Liczby losowe generowane przez programy komputerowe nazywamy liczbami pseudolosowymi.

Dysponując liczbami losowymi $u_i$ należącymi do rozkładu równomiernego $U(0,1)$ możemy łatwo uzyskać liczby losowe $x_i$ należące do rozkładu równomiernego $U(a,b)$. Wystarczy skorzystać ze wzoru:

$x_i=a+(b-a)u_i$

Przykład histogramu rozkładu równomiernego $U(0,1)$ Czytelnik znajdzie w arkuszu kalkulacyjnym: Równomierny

Dystrybuanta odwrotna

Dystrybuanta odwrotna $F^{-1}(u)$ przekształca zmienną losową $U(0,1)$ o rozkładzie równomiernym w zmienną losową $X$ o rozkładzie $f(x)$, któremu odpowiada dystrybuanta $F(x)$:

$X=F^{-1}(U)$

Punktowi $u\in [0,1]$ zostaje przypisany punkt $x=F^{-1}(u)\in [a,b]$. Granice przedziału $[a,b]$ mogą być w ogólności niewłaściwe.

Przykładowo dystrybuanta odwrotna zmiennej losowej o rozkładzie $f(x)=e^{-x}$, gdzie $x\in [0,+\mathrm{\infty })$ jest równa $F^{-1}(x)=-\log (1-x)$. Zatem wylosowana liczba $u_1=0.25\in [0.2,0.3)$ z rozkładu $U[0,1]$ przejdzie w punkt

$x_1=-\log (1-u_1)=0.125\in [0.1,0.2)$

z rozkładu wykładniczego $f(x)=e^{-x}$. Podobnie liczba $u_2=0.95\in [0.9,0.1)$ przejdzie w punkt

$x_2=-\log (1-u_2)=1.301\in [1.3,1.4)$

Czyli liczby te trafią do innych podprzedziałów, będą zliczane w innych miejscach i utworzą inny histogram. Tak jak to pokazano na rysunku:

Przykłady histogramów rozkładu równomiernego $U(0,1)$ i wykładniczego $f(x)=e^{-x}$ (wygenerowanego z rozkładu równomiernego) Czytelnik znajdzie w arkuszach kalkulacyjnych:

Równomierny
Wykładniczy

Przedstawimy teraz kilka przykładów zastosowania tego faktu.

Rozkład jednomianowy $f(x)=(n+1)x^n$

Rozważmy funkcję gęstości prawdopodobieństwa postaci:

$f(x)=\{\begin{array}{lll}0& \text{dla}& x<0\\ (n+1)x^n& \text{dla}& x\in [0,1]\\ 0& \text{dla}& x>1\end{array}$

Łatwo sprawdzamy, że

${\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}f(t)dt={\int }_0^{1}f(t)dt=(n+1){\int }_0^1{t}^{n}dt=t^{n+1}{|}_0^1=1$

Dystrybuanta

$F(x)=\{\begin{array}{lll}0& \text{dla}& x<0\\ x^{n+1}& \text{dla}& x\in [0,1]\\ 1& \text{dla}& x>1\end{array}$

Dystrybuanta odwrotna

$F^{-1}(x)=\sqrt[n+1]{x}$

Jeżeli $u_i\in U(0,1)$, to liczby $x_i=\sqrt[n+1]{u_i}\in [0,1]$ będą należały do rozkładu $f(x)=(n+1)x^n$ określonego na odcinku $[0,1]$.

Krzywa opisująca histogram

Załóżmy, że:

• wygenerowaliśmy $N$ liczb losowych $u_i\in U(0,1)$

• obliczyliśmy wartości $x_i=\sqrt[n+1]{u_i}$

• podzieliliśmy przedział zmienności liczb $x_i$ (w naszym przypadku $[0,1]$) na podprzedziały każdy o ustalonej szerokości $\mathrm{\Delta }$

• pogrupowaliśmy $x_i$ w poszczególnych podprzedziałach i wyznaczyliśmy ilość $g(k)$ liczb $x_i$ w $k$-tym podprzedziale

Jakiej zależności $g(k)$ należy się spodziewać? Prawdopodobieństwo tego, że zmienna losowa przyjmie wartość z $k$-tego przedziału o szerokości $\mathrm{\Delta }$ wyraża się wzorem:

$P((k-1)\mathrm{\Delta }⩽x_i⩽k\mathrm{\Delta })={\int }_{(k-1)\mathrm{\Delta }}^{k\mathrm{\Delta }}f(t)dt\approx f(k\mathrm{\Delta })\cdot \mathrm{\Delta }$

Przybliżenie jest tym lepsze, im mniejsza jest szerokość przedziałów $\mathrm{\Delta }$. Należy zatem oczekiwać zależności:

$g(k)=N\cdot f(k\mathrm{\Delta })\cdot \mathrm{\Delta }=N\mathrm{\Delta }\cdot f(k\mathrm{\Delta })$

Dla rozkładu jednomianowego na odcinku $[0,1]$ otrzymamy:

$g(k)=N\cdot f(k\mathrm{\Delta })\cdot \mathrm{\Delta }=N\cdot (n+1)\cdot (k\mathrm{\Delta }{)}^n\cdot \mathrm{\Delta }=(n+1)N\cdot {\mathrm{\Delta }}^{n+1}\cdot k^n$

Przykład dla $n=1$ i $n=2$

Niech $N={10}^5$, $\mathrm{\Delta }=0.01$

dla $n=1$: $g(k)=20k$

dla $n=2$: $g(k)=0.3k^2$

Przykłady histogramów rozkładu jednomianowego $f(x)=(n+1)x^n$, dla $n=1,2$ Czytelnik znajdzie w arkuszach kalkulacyjnych:

Jednomianowy (n = 1)
Jednomianowy (n = 2)

Rozkład postaci $f(x)=\frac{3}{2}\sqrt{x}$

Rozważmy następującą funkcję gęstości prawdopodobieństwa:

$f(x)=\{\begin{array}{lll}0& \text{dla}& x<0\\ \frac{3}{2}\sqrt{x}& \text{dla}& x\in [0,1]\\ 0& \text{dla}& x>1\end{array}$

Czytelnik łatwo sprawdzi, że:

${\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}f(t)dt=1$

Dystrybuanta

$F(x)=\{\begin{array}{lll}0& \text{dla}& x<0\\ x^{3/2}& \text{dla}& x\in [0,1]\\ 1& \text{dla}& x>1\end{array}$

Dystrybuanta odwrotna

$F^{-1}(x)=x^{2/3}$

Jeżeli $u_i\in U(0,1)$, to liczby $x_i=(u_i{)}^{2/3}$ będą należały do rozkładu $f(x)=\frac{3}{2}\sqrt{x}$ określonego na odcinku $[0,1]$.

Krzywa opisująca histogram

$g(k)=N\cdot f(k\mathrm{\Delta })\cdot \mathrm{\Delta }=N\cdot \frac{3}{2}\cdot \sqrt{k\mathrm{\Delta }}\cdot \mathrm{\Delta }=\frac{3}{2}\cdot N{\mathrm{\Delta }}^{3/2}\cdot \sqrt{k}$

Dla $N={10}^5$, $\mathrm{\Delta }=0.01$ otrzymujemy: $g(k)=150\sqrt{k}$

Przykład histogramu rozkładu postaci $f(x)=\frac{3}{2}\sqrt{x}$ Czytelnik znajdzie w arkuszu kalkulacyjnym: Pierwiastek

Rozkład postaci $f(x)=\frac{1}{2\sqrt{ax}}$

Dla $a>0$ określamy następującą funkcję gęstości prawdopodobieństwa:

$f(x)=\{\begin{array}{lll}0& \text{dla}& x⩽0\\ \frac{1}{2\sqrt{ax}}& \text{dla}& x\in (0,a]\\ 0& \text{dla}& x>a\end{array}$

Czytelnik łatwo sprawdzi, że:

${\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}f(t)dt=1$

Dystrybuanta

$F(x)=\{\begin{array}{lll}0& \text{dla}& x<0\\ \sqrt{\frac{x}{a}}& \text{dla}& x\in [0,a]\\ 1& \text{dla}& x>a\end{array}$

Dystrybuanta odwrotna

$F^{-1}(x)=a{x}^2$

Jeżeli $u_i\in U(0,1)$, to liczby $x_i=a{u}_i^2\in [0,a]$ będą należały do rozkładu $f(x)=\frac{1}{2\sqrt{ax}}$ określonego na odcinku $(0,a]$.

Krzywa opisująca histogram

$g(k)=N\cdot f(k\mathrm{\Delta })\cdot \mathrm{\Delta }=N\cdot \frac{1}{2\sqrt{ak\mathrm{\Delta }}}\cdot \mathrm{\Delta }=\frac{N\sqrt{\mathrm{\Delta }}}{2\sqrt{a}}\cdot \frac{1}{\sqrt{k}}$

Dla $a=25$, $N={10}^5$, $\mathrm{\Delta }=0.25$ otrzymujemy: $g(k)=\frac{5000}{\sqrt{k}}$

Przykład histogramu rozkładu postaci $f(x)=\frac{1}{10\sqrt{x}}$ określonego na odcinku $(0,25]$ Czytelnik znajdzie w arkuszu kalkulacyjnym: Odwrotność pierwiastka

Rozkład wykładniczy $f(x)=\lambda e^{-\lambda x}$

Dla rozkładu wykładniczego funkcja gęstości prawdopodobieństwa jest określona następująco:

$f(x)=\{\begin{array}{lll}0& \text{dla}& x<0\\ \lambda e^{-\lambda x}& \text{dla}& x⩾0\end{array}$

gdzie $\lambda >0$.

Dystrybuanta

$F(x)=\{\begin{array}{lll}0& \text{dla}& x<0\\ 1-e^{-\lambda x}& \text{dla}& x⩾0\end{array}$

Dystrybuanta odwrotna

$F^{-1}(x)=-\frac{1}{\lambda }\cdot \log (1-x)$

Jeżeli $u_i\in U(0,1)$, to liczby $x_i=-\frac{1}{\lambda }\cdot \log (1-u_i)$ będą należały do rozkładu wykładniczego $f(x)=\lambda e^{-\lambda x}$ określonego na półprostej $[0,+\mathrm{\infty })$.

Krzywa opisująca histogram

$g(k)=N\cdot \lambda e^{-\lambda \cdot k\mathrm{\Delta }}\cdot \mathrm{\Delta }=N\mathrm{\Delta }\cdot \exp (-(\lambda \mathrm{\Delta })\cdot k)$

Dla $\lambda =1$ oraz $N={10}^5$, $\mathrm{\Delta }=0.1$ otrzymujemy: $g(k)=10000\cdot e^{-0.1\cdot k}$

Przykład histogramu rozkładu wykładniczego $f(x)=e^{-x}$ Czytelnik znajdzie w arkuszu kalkulacyjnym: Wykładniczy

Rozkład normalny $N(\mu ,{\sigma }^2)$

Rozkładem normalnym nazywamy rozkład, dla którego funkcja gęstości prawdopodobieństwa ma postać:

$f(x)=\frac{1}{\sigma \sqrt{2\pi }}\cdot \exp (-\frac{(x-\mu {)}^2}{2{\sigma }^2})$

gdzie $\mu \in \mathbb{R}$ i $\sigma >0$.

Dystrybuanta

$F(x)={\int }_{-\mathrm{\infty }}^{x}f(t)dt=\frac{1}{\sqrt{\pi }}{\int }_{-\mathrm{\infty }}^{\frac{x-\mu }{\sigma \sqrt{2}}}{e}^{-u^2}du=\frac{1}{\sqrt{\pi }}{\int }_{-\mathrm{\infty }}^0{e}^{-u^2}du+\frac{1}{\sqrt{\pi }}{\int }_0^{\frac{x-\mu }{\sigma \sqrt{2}}}{e}^{-u^2}du$

Ponieważ

$\frac{1}{\sqrt{\pi }}{\int }_{-\mathrm{\infty }}^0{e}^{-u^2}du=\frac{1}{2}$

to

$F(x)=\frac{1}{2}(1+\text{erf}\left(\frac{x-\mu }{\sigma \sqrt{2}}\right))$

Funkcję $\text{erf}(x)$ nazywamy funkcją błędu Gaussa i jest to funkcja nieelementarna:

$\text{erf}(x)=\frac{2}{\sqrt{\pi }}{\int }_0^x{e}^{-t^2}dt$

Łatwo można pokazać, że $\text{erf}(x)$ jest funkcją nieparzystą:

$\text{erf}(-x)=-\text{erf}(x)$

W arkuszu LibreOffice $\text{erf}(x)$ jest dostępna pod nazwą FUNKCJA.BŁ.DOKŁ(). Dostępna jest też komplementarna funkcja błędu $\text{erfc}(x)=1-\text{erf}(x)$ pod nazwą KOMP.FUNKCJA.BŁ.DOKŁ(). Funkcja odwrotna funkcji $\text{erf}(x)$ również nie jest elementarna. Dlatego uzyskanie liczb $x_i\in N(\mu ,{\sigma }^2)$ na bazie liczb $u_i\in U(0,1)$ wymaga nieco innego podejścia.

Metoda Boxa - Mullera

Zamiast jednej zmiennej losowej o standardowym rozkładzie normalnym $N(0,1)$, rozważmy dwie niezależne zmienne losowe o takim rozkładzie. Funkcja gęstości prawdopodobieństwa dla takiej pary niezależnych zmiennych losowych będzie iloczynem gęstości prawdopodobieństwa tych funkcji:

$f(x,y)=f(x)f(y)=\frac{1}{\sqrt{2\pi }}{e}^{-x^2/2}\cdot \frac{1}{\sqrt{2\pi }}{e}^{-x^2/2}=\frac{1}{2\pi }{e}^{-(x^2+y^2)/2}$

Przechodząc do współrzędnych biegunowych

$\{\begin{array}{l}x=r\cdot \cos (\phi )\\ y=r\cdot \sin (\phi )\end{array}$

gdzie $r⩾0$ i $\phi \in [0,2\pi )$, otrzymujemy:

$f(r,\phi )=\frac{1}{2\pi }{e}^{-r^2/2}$

Funkcję gęstości prawdopodobieństwa $f(r,\phi )$ możemy zapisać w postaci iloczynu:

$f(r,\phi )=g(r)h(\phi )$

Widzimy, że $h(\phi )=\frac{1}{2\pi }=\text{const}$ jest unormowaną funkcją gęstości prawdopodobieństwa zmiennej $\phi$. Oznacza to, że rozkład $\phi$ jest rozkładem równomiernym $U(0,2\pi )$, a liczba $\phi$ może być zapisana w postaci $\phi =2\pi \cdot u$, gdzie $u$ jest liczbą z równomiernego rozkładu $U(0,1)$.

Iloczyn dystrybuant $G(r)H(\phi )$ jest określony całką we współrzędnych biegunowych:

$G(\hat{r})H(\hat{\phi })={\int }_0^{\hat{r}}{\int }_0^{\hat{\phi }}g(r)h(\phi )rdrd\phi ={\int }_0^{\hat{r}}g(r)rdr{\int }_0^{\hat{\phi }}h(\phi )d\phi$

Zatem dystrybuanta $G(\hat{r})$ jest równa:

$G(\hat{r})={\int }_0^{\hat{r}}g(r)rdr={\int }_0^{\hat{r}}{e}^{-r^2/2}rdr=-e^{-r^2/2}{|}_0^{\hat{r}}=1-e^{-{\hat{r}}^2/2}$

Całkę nieoznaczoną $\int e^{-r^2/2}rdr=-e^{-r^2/2}$ wyliczamy dokonując podstawienia $t=-r^2/2$.

Wracając do zmiennej $r$, mamy:

$G(r)=1-e^{-r^2/2}$

Łatwo znajdujemy dystrybuantę odwrotną:

$G^{-1}(r)=\sqrt{-2\log (1-r)}$

Jeżeli $v_i\in U(0,1)$, to liczby $\sqrt{-2\log (1-v_i)}$ będą należały do rozkładu prawdopodobieństwa odpowiadającego funkcji $g(r)$.

Zatem parze liczb $u,v\in U(0,1)$ odpowiadają liczby $\phi ,r$

$\{\begin{array}{l}\phi =2\pi \cdot u\\ r=\sqrt{-2\log (1-v)}\end{array}$

z rozkładów opisywanych funkcjami gęstości $h(\phi )$ i $g(r)$, a tym liczbom odpowiada para liczb $x,y$

$\{\begin{array}{l}x=\sqrt{-2\log (1-v)}\cdot \cos (2\pi \cdot u)\\ y=\sqrt{-2\log (1-v)}\cdot \sin (2\pi \cdot u)\end{array}$

które należą do standardowych rozkładów normalnych $N(0,1)$.

Wnioski

Jeżeli $u,v\in U(0,1)$, to liczby

$\{\begin{array}{l}x=\sqrt{-2\log (1-u)}\cdot \cos (2\pi \cdot v)\\ y=\sqrt{-2\log (1-u)}\cdot \sin (2\pi \cdot v)\end{array}$

będą należały do standardowego rozkładu normalnego $N(0,1)$.

Uogólniając postępowanie z poprzedniego punktu, można łatwo pokazać, że jeżeli $u,v\in U(0,1)$, to liczby

$\{\begin{array}{l}x=\mu +\sigma \sqrt{-2\log (1-u)}\cdot \cos (2\pi \cdot v)\\ y=\mu +\sigma \sqrt{-2\log (1-u)}\cdot \sin (2\pi \cdot v)\end{array}$

będą należały do rozkładu normalnego $N(\mu ,{\sigma }^2)$.

Przykład histogramu rozkładu normalnego $f(x)=\frac{1}{\sqrt{2\pi }}{e}^{-x^2}$ Czytelnik znajdzie w arkuszu kalkulacyjnym: Normalny