Liczba e ukryta w sumie rozkładów jednostajnych

Rozkład jednostajny na odcinku (0,1), chyba najprostszy z możliwych rozkładów ciągłych, z pozoru niezbyt interesujący, a jednak 🙂 Dziś ciekawostka wiążąca rozkład sumy rozkładów jednostajnych z liczbą Eulera e.

Uniform Sum Distribution

Rozkład jednostajny ciągły na odcinku (a,b)

Rozkład jednostajny ciągły na odcinku (a,b) jest opisany poniższą funkcją gęstości.

f(x)=\begin{cases}\frac{1}{b-a}&&\text{dla }a\leq x\leq b\\0&&\text{w p.p.}\end{cases}

Pisząc X\sim U(a,b) oznaczamy, że zmienna losowa X ma rozkład jednostajny ciągły na odcinku (a,b). Jest to rozkład ciągły, zatem przyjęcie wartości 0 lub \frac{1}{b-a} w punktach x=a i x=b jest umowne i nie ma zwykle wpływu na własności i rozważania.

Rozkład jednostajny - gęstość

Rozkład jednostajny przypisuje prawdopodobieństwo zdarzenia na bazie jego "rozmiaru", a dokładniej rozmiaru tej części, która znajduje się w odcinku (a,b).

Jeśli X\sim U(a,b) to:

  • \rm E X=\frac{a+b}{2}
  • \rm Mediana(X)=\frac{a+b}{2}
  • \rm Var(X)=\frac{(b-a)^2}{12}

Rozkład Irwina-Halla - czyli suma rozkładów jednostajnych ciągłych

Dla n niezależnych zmiennych losowych pochodzących z rozkładu jednostajnego ciągłego na odcinku (0,1)

X_1,X_2,\ldots,X_n\sim U(0,1)

można zdefiniować zmienną losową będącą ich sumą

X=X_1+X_2+\ldots+X_n=\displaystyle\sum_{i=1}^nX_i

O zmiennej X mówimy, że pochodzi z rozkładu Irwina-Halla, co zapisujemy

X\sim IH(n)

Funkcja gęstości rozkładu Irwina-Halla dana jest poniższym wzorem (wyprowadzenie pomijam, gdyż nie jest to główny temat wpisu).

f(x)=\frac{1}{(n-1)!}\displaystyle\sum_{k=0}^{\lfloor x\rfloor}(-1)^k{n\choose k}(x-k)^{n-1}

Rozkład Irwina-Halla

Jeśli X\sim IH(n) to:

  • \rm E X=\frac{n}{2}
  • \rm Mediana(X)=\frac{n}{2}
  • \rm Var(X)=\frac{n}{12}

Ciekawostka: rozkład sumy rozkładów jednostajnych przypomina "kształtem" rozkład normalny. W praktyce możliwe jest "przybliżone" wygenerowanie liczb z rozkładu normalnego korzystając z niezbyt "długiej" sumy liczb z rozkładu jednostajnego, odpowiednio przesuwając i skalując wynik.

Jak "standaryzować" funkcję gęstości?

  • f(x) - gęstość rozkładu - znana
  • EX=\mu - znane
  • Var(X)=\sigma^2 - znane

Poszukujemy nowej gęstości rozkładu prawdopodobieństwa

g(x)=af(bx+c)

Zaczynamy od warunku na gęstość - tzn. "masa prawdopodobieństwa" to 1.

\displaystyle\int_\mathbb{R}g(x)dx=1

\displaystyle\int_\mathbb{R}g(x)dx=\displaystyle\int_\mathbb{R}af(bx+c)dx=

=\begin{bmatrix}t=bx+c\\\frac{dt}{b}=dx\end{bmatrix}=

=\displaystyle\int_\mathbb{R}af(t)\frac{dt}{b}=\frac{a}{b}\displaystyle\int_\mathbb{R}f(t)dt=\frac{a}{b}\cdot 1

\frac{a}{b}=1 zatem a=b

Zapisujemy

g(x)=af(ax+c)

Stosujemy warunek na wartość oczekiwaną równą "0"

0=\displaystyle\int_\mathbb{R}xg(x)dx=\displaystyle\int_\mathbb{R}xaf(ax+c)dx=

=\begin{bmatrix}t=ax+c\\\frac{dt}{a}=dx\\x=\frac{t-c}{a}\end{bmatrix}=

=\displaystyle\int_\mathbb{R}\frac{t-c}{a}af(t)\frac{dt}{a}=\frac{1}{a}\displaystyle\int_\mathbb{R}(t-c)f(t)dt=

=\frac{1}{a}\Bigg[\displaystyle\int_\mathbb{R}tf(t)dt-c\displaystyle\int_\mathbb{R}f(t)dt\Bigg]=\frac{1}{a}(\mu-c)

\frac{\mu-c}{a}=0

\mu-c=0 zatem c=\mu

Zapisujemy

g(x)=af(ax+\mu)

Stosujemy warunek na wariancję równą "1"

1=\displaystyle\int_\mathbb{R}x^2g(x)dx-\Bigg(~~\underbrace{\displaystyle\int_\mathbb{R}xg(x)dx}_{0}~~\Bigg)^2=\displaystyle\int_\mathbb{R}x^2af(ax+\mu)dx

=\begin{bmatrix}t=ax+\mu\\\frac{dt}{a}=dx\\x=\frac{t-\mu}{a}\end{bmatrix}=

=\displaystyle\int_\mathbb{R}\Bigg(\frac{t-\mu}{a}\Bigg)^2af(t)\frac{dt}{a}=\frac{1}{a^2}\displaystyle\int_\mathbb{R}\Big(t^2-2t\mu+\mu^2\Big)f(t)dt=

=\frac{1}{a^2}\Bigg[\displaystyle\int_\mathbb{R}t^2f(t)dt-2\mu\displaystyle\int_\mathbb{R}tf(t)dt+\mu^2\displaystyle\int_\mathbb{R}f(t)dt\Bigg]=...

ale

\sigma^2=\displaystyle\int_\mathbb{R}t^2f(t)dt-\Bigg(\displaystyle\int_\mathbb{R}tf(t)dt\Bigg)^2=\displaystyle\int_\mathbb{R}t^2f(t)dt-\mu^2

więc

\displaystyle\int_\mathbb{R}t^2f(t)dt=\sigma^2+\mu^2

...=\frac{1}{a^2}\Big(\sigma^2+\mu^2-2\mu^2+\mu^2\Big)=\frac{\sigma^2}{a^2}

\frac{\sigma^2}{a^2}=1

a^2=\sigma^2

a=\sigma

Funkcja gęstości jest funkcją dodatnią, dlatego odrzucamy a=-\sigma.

Ostatecznie

g(x)=\sigma f(\sigma x+\mu)

Podsumowując, jeśli zmienna losowa X pochodzi z rozkładu opisanego funkcją gęstości f(x), o wartości oczekiwanej EX=\mu oraz skończonej wariancji Var(X)=\sigma^2, to standaryzowana zmienna losowa

\frac{X-\mu}{\sigma}

pochodzi z rozkładu opisanego gęstością

\sigma f(\sigma x+\mu)

Jakość przybliżenia rozkładu normalnego na bazie sumy rozkładów jednostajnie ciągłych

Porównany dwie funkcje gęstości

1. Gęstość rozkładu N(0,1) tzw. standardowego rozkładu normalnego

\frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}}

2. Standaryzowana gęstość rozkładu Irwina-Halla

\sqrt{\frac{n}{12}}\frac{1}{(n-1)!}\displaystyle\sum_{k=0}^{\bigg\lfloor x\sqrt{\frac{n}{12}}+\frac{n}{2}\bigg\rfloor}(-1)^k{n\choose k}\Bigg(x\sqrt{\frac{n}{12}}+\frac{n}{2}-k\Bigg)^{n-1}

Poniżej animacja obrazująca bardzo szybką zbieżność gęstości do N(0,1)

Rozkład normalny jako suma rozkładów jednostajnych

Gdzie ten Euler?

Korzystając z centralnego twierdzenia granicznego, niemal na bazie każdego rozkładu, jesteśmy w stanie wygenerować liczby losowe przybliżające rozkład normalny - zatem, w tym sensie, rozkład Irwina-Halla niczym się nie wyróżnia. Jest w nim jednak coś specjalnego - szczególne powiązanie z liczbą Eulera e. W dalszej części tekstu będą rozważał minimalną liczbę losowań z rozkładu U(0,1) niezbędnych do przekroczenia pewnej ustalonej wartości - np. losujemy, sumujemy, kończymy po otrzymaniu 1 lub więcej.

Minimalna liczba losowań z U(0,1) potrzebna do przekroczenia liczby 1 (przekroczenie osiągane sumując wyniki)

Przez N_1 oznaczamy taką zmienną losową powstającą na bazie n niezależnych zmiennych losowych z rozkładu jednostajnego ciągłego U(0,1), że

N_1=\min\Big\{n~:~X_1+X_2+\ldots+X_n\geq 1\Big\}

X_i\sim U(0,1)

Jeśli N_1=n to X_1+X_2+\ldots+X_n\geq 1 oraz X_1+X_2+\ldots+X_{n-1}<1.

Analogicznie tworzymy zmienną losową N_a

N_a=\min\Big\{n~:~X_1+X_2+\ldots+X_n\geq a\Big\}

X_i\sim U(0,1)

Jeśli N_a=n to X_1+X_2+\ldots+X_n\geq a oraz X_1+X_2+\ldots+X_{n-1}<a.

Wartość oczekiwania zmiennej losowej N_1?

Spostrzeżenie #1: N_1 jest zmienną losową dyskretną o niezbyt "przyjemnej" w analizie funkcji masy prawdopodobieństwa.

pmf_{N_1}(n) - oznacza prawdopodobieństwo, że dokładnie n losowań wystarczy do przekroczenia jedności, natomiast n-1 losowań to zbyt mało.

Spostrzeżenie #2: dystrybuanta rozkładu, z którego pochodzi N_1, jest łatwiejsza w analizie, gdyż podaje prawdopodobieństwo, że n losowań wystarczy (bez żadnych dodatkowych warunków).

F_{N_1}(n)=P\Big(N_1\leq n\Big)

Spostrzeżenie #3: zmienna losowa N_1 przyjmuje wartości nieujemne, dlatego jej wartość oczekiwaną można wyrazić za pomocą dystrybuanty (sumując pole powierzchni ponad dystrybuantą).

EN_1=\displaystyle\sum_{n=0}^\infty\Big(1-F_{N_1}(n)\Big)

Spostrzeżenie #4: P\Big(N_1=0\Big)=0 oraz P\Big(N_1=1\Big)=0 zatem F_{N_1}(0)=0 oraz F_{N_1}(1)=0, wtedy

EN_1=1+1+\displaystyle\sum_{n=2}^\infty\Big(1-F_{N_1}(n)\Big)

Spostrzeżenie #5: jeśli F_{N_1}(n) oznacza prawdopodobieństwo, że n losowań wystarczy, to 1-F_{N_1}(n) oznacza prawdopodobieństwo, że n losowań to zbyt mało.

Spostrzeżenie #6: w przestrzeni \mathbb{R}^n zdarzenie, że n losowań to zbyt mało, reprezentowane jest przez zbiór

\bigg\{(x_1,x_2,\ldots,x_n)\in\mathbb{R}^n~:~0\leq x_1+x_2+\ldots+x_n<1\bigg\}

Bez wpływu na wartość prawdopodobieństwa (dodajemy jedynie "górną powierzchnię") możemy zapisać

\bigg\{(x_1,x_2,\ldots,x_n)\in\mathbb{R}^n~:~0\leq x_1+x_2+\ldots+x_n\leq 1\bigg\}

Spostrzeżenie #7: graficzna reprezentacja ww. zbiorów dla przestrzeni 1, 2 i 3 wymiarowej ujawnia prawidłowość.

1-sympleks

2-sympleks

3-sympleks

Spostrzeżenie #8: rozważany zbiór

\bigg\{(x_1,x_2,\ldots,x_n)\in\mathbb{R}^n~:~0\leq x_1+x_2+\ldots+x_n\leq 1\bigg\}

to n-sympleks.

N-sympleks i jego objętość

1, 2, 3 - sympleks
Sympleks w matematyce to uogólnienie pojęcia trójkąta na przestrzenie wielowymiarowe. Mówimy, że n-sympleks powstaje z n+1 wierzchołków, będąc obiektem osadzonym w przestrzeni k-wymiarowej. Sam sympleks to "otoczka" wypukła jego wierzchołków, jest to więc "powierzchnia", jeśli określenie "powierzchni" ma jakikolwiek sens dla przestrzeni wielowymiarowych.

  • Niech \lambda będzie miarą Lebesgue’a w \mathbb{R}^n - miara Lebesgue’a jest wielowymiarowym uogólnieniem pojęcia długości, pola powierzchni, objętości.
  • Niech V_n(a) oznacza objętość n-sympleksu przeskalowanego współczynnikiem skali podobieństwa a - wtedy V_n(1) to objętość "standardowego" n-sympleksu.

Możemy zapisać

V_n(1)=\lambda\Bigg(\Big\{0\leq x_1+x_2+\ldots+x_n\leq 1\Big\}\Bigg)

V_n(a)=\lambda\Bigg(\Big\{0\leq x_1+x_2+\ldots+x_n\leq a\Big\}\Bigg)

Zgodnie z założeniem, że a jest współczynnikiem skali podobieństwa w przestrzeni n-wymiarowej

V_n(a)=a^nV_n(1)

V_n(1)=\displaystyle\int_{x_1+x_2+\ldots+x_n\leq 1}dx_1dx_2\ldots dx_n

V_n(a)=\displaystyle\int_{x_1+x_2+\ldots+x_n\leq a}dx_1dx_2\ldots dx_n

Zauważmy, że x_1+x_2+\ldots+x_n\leq 1 wtedy i tylko wtedy gdy x_n\leq 1 oraz x_1+x_2+\ldots+x_{n-1}\leq 1-x_n, więc

V_n(1)=\displaystyle\int_{x_n\leq 1}\Bigg(~~\underbrace{\displaystyle\int_{x_1+x_2+\ldots+x_{n-1}\leq 1-x_n}dx_1dx_2\ldots dx_{n-1}}_{V_{n-1}(1-x_n)}~~\Bigg)dx_n=

=\displaystyle\int_{x_n\leq 1}V_{n-1}(1-x_n)dx_n=\displaystyle\int_{x_n\leq 1}(1-x_n)^{n-1}V_{n-1}(1)dx_n=

=V_{n-1}(1)\displaystyle\int_{0}^1(1-x_n)^{n-1}dx_n=

=\begin{bmatrix}t=1-x_n&&t_0=1-0=1\\-dt=dx_n&&t_1=1-1=0\end{bmatrix}=

=V_{n-1}(1)\Bigg(-\displaystyle\int_{1}^0t^{n-1}dt\Bigg)=

=V_{n-1}(1)\displaystyle\int_{0}^1t^{n-1}dt=V_{n-1}(1)\Bigg[\frac{t^n}{n}\Bigg]_0^1=\frac{1}{n}V_{n-1}(1)

Finalnie otrzymaliśmy wzór rekurencyjny

V_n(1)=\frac{1}{n}V_{n-1}(1)

V_1(1)=1

Podstawiając kolejne wartości

V_2(1)=\frac{1}{2}\cdot 1=\frac{1}{1\cdot 2}

V_3(1)=\frac{1}{3}\cdot\frac{1}{1\cdot 2}=\frac{1}{1\cdot 2\cdot 3}

ujawnia się wzór na objętość n-sympleksu.

V_n(1)=\frac{1}{n!}

V_n(a)=\frac{a^n}{n!}

Powrót do wartości oczekiwanej EN_1

1-F_{N_1}(n)=V_n(1)=\frac{1}{n!}

EN_1=1+1+\displaystyle\sum_{n=2}^\infty\Big(1-F_{N_1}(n)\Big)=\frac{1}{0!}+\frac{1}{1!}+\displaystyle\sum_{n=2}^\infty\frac{1}{n!}

Ostatecznie

EN_1=\displaystyle\sum_{n=0}^\infty\frac{1}{n!}=e

Wartości oczekiwane EN_1EN_2EN_3EN_4EN_5

  • EN_1=e
  • EN_2=e^2-e
  • EN_3=e^3-2e^2+\frac{e}{2}
  • EN_4=e^4-3e^3+2e^2-\frac{e}{6}
  • EN_5=e^5-4e^4+\frac{9e^3}{2}-\frac{4e^2}{3}+\frac{e}{24}

Po szczegóły zapraszam do artykułu "Uniform Sum Distribution".

Weryfikacja na bazie Monte Carlo

Eksperyment wykonany przy wykorzystaniu MathParser.org-mXparser.


/*
* Rekurencja w poszukiwaniu minimalnego 'n'
* a - liczba do przekroczenia
* s - tymczasowa suma
* k - numer losowania
*/
Function Nrec = new Function("Nrec(a, s, k) = if( s >= a, k, Nrec( a, s + [Uni], k+1 ) )");
/*
* Definicja zmiennych losowych
*/
Argument N1 = new Argument("N1 = Nrec(1,0,0)", Nrec);
Argument N2 = new Argument("N2 = Nrec(2,0,0)", Nrec);
Argument N3 = new Argument("N3 = Nrec(3,0,0)", Nrec);
Argument N4 = new Argument("N4 = Nrec(4,0,0)", Nrec);
Argument N5 = new Argument("N5 = Nrec(5,0,0)", Nrec);
/*
* Estymacja wartości oczekiwanych - 1000000 operacji
*/
Expression EN1 = new Expression("avg( i, 1, 1000000, N1 )", N1);
Expression EN2 = new Expression("avg( i, 1, 1000000, N2 )", N2);
Expression EN3 = new Expression("avg( i, 1, 1000000, N3 )", N3);
Expression EN4 = new Expression("avg( i, 1, 1000000, N4 )", N4);
Expression EN5 = new Expression("avg( i, 1, 1000000, N5 )", N5);
/*
* Teoretyczne wartości oczekiwane
*/
Expression eN1 = new Expression("e");
Expression eN2 = new Expression("e^2 - e");
Expression eN3 = new Expression("e^3 - 2*e^2 + e/2");
Expression eN4 = new Expression("e^4 - 3*e^3 + 2*e^2 - e/6");
Expression eN5 = new Expression("e^5 - 4*e^4 + 9/2 * e^3 - 4/3 * e^2 + e/24");
/*
* Obliczenie i wyswietlenie wyników
*/
mXparser.consolePrintln("EN1 - estymowane: " + EN1.getExpressionString() + " = " + EN1.calculate());
mXparser.consolePrintln("EN1 - teoretyczne: " + eN1.getExpressionString() + " = " + eN1.calculate());
mXparser.consolePrintln("------------------------------");
mXparser.consolePrintln("EN2 - estymowane: " + EN2.getExpressionString() + " = " + EN2.calculate());
mXparser.consolePrintln("EN2 - teoretyczne: " + eN2.getExpressionString() + " = " + eN2.calculate());
mXparser.consolePrintln("------------------------------");
mXparser.consolePrintln("EN3 - estymowane: " + EN3.getExpressionString() + " = " + EN3.calculate());
mXparser.consolePrintln("EN3 - teoretyczne: " + eN3.getExpressionString() + " = " + eN3.calculate());
mXparser.consolePrintln("------------------------------");
mXparser.consolePrintln("EN4 - estymowane: " + EN4.getExpressionString() + " = " + EN4.calculate());
mXparser.consolePrintln("EN4 - teoretyczne: " + eN4.getExpressionString() + " = " + eN4.calculate());
mXparser.consolePrintln("------------------------------");
mXparser.consolePrintln("EN5 - estymowane: " + EN5.getExpressionString() + " = " + EN5.calculate());
mXparser.consolePrintln("EN5 - teoretyczne: " + eN5.getExpressionString() + " = " + eN5.calculate());
mXparser.consolePrintln("------------------------------");

Wynik


[mXparser-v.4.1.1] EN1 - estymowane: avg( i, 1, 1000000, N1 ) = 2.71813
[mXparser-v.4.1.1] EN1 - teoretyczne: e = 2.718281828459045
[mXparser-v.4.1.1] ------------------------------
[mXparser-v.4.1.1] EN2 - estymowane: avg( i, 1, 1000000, N2 ) = 4.668181
[mXparser-v.4.1.1] EN2 - teoretyczne: e^2 - e = 4.670774270471606
[mXparser-v.4.1.1] ------------------------------
[mXparser-v.4.1.1] EN3 - estymowane: avg( i, 1, 1000000, N3 ) = 6.666433
[mXparser-v.4.1.1] EN3 - teoretyczne: e^3 - 2*e^2 + e/2 = 6.666565639555883
[mXparser-v.4.1.1] ------------------------------
[mXparser-v.4.1.1] EN4 - estymowane: avg( i, 1, 1000000, N4 ) = 8.664254
[mXparser-v.4.1.1] EN4 - teoretyczne: e^4 - 3*e^3 + 2*e^2 - e/6 = 8.66660449003271
[mXparser-v.4.1.1] ------------------------------
[mXparser-v.4.1.1] EN5 - estymowane: avg( i, 1, 1000000, N5 ) = 10.667831
[mXparser-v.4.1.1] EN5 - teoretyczne: e^5 - 4*e^4 + 9/2 * e^3 - 4/3 * e^2 + e/24 = 10.66666206862245

Wszystko się zgadza! Niesamowite powiązanie!

Pozdrowienia,

Mariusz Gromada

Views All Time
Views All Time
609
Views Today
Views Today
1

Dodaj komentarz

Twój adres email nie zostanie opublikowany. Pola, których wypełnienie jest wymagane, są oznaczone symbolem *