Rachunek prawdopodobieństwa i statystyka/Wykład 14: Komputerowe metody statystyki

From Studia Informatyczne

Omówimy, bardzo krótko, trzy zagadnienia związane z możliwością zastosowania metod informatycznych w statystyce: generowanie liczb o charakterze losowym, metodę bootstrap oraz jądrową estymację gęstości.

Na ostatnim wykładzie poznamy pewne metody, które zostały opracowane w drugiej połowie XX wieku, a ich przewodnią ideą jest wykorzystanie dużej mocy obliczeniowej komputera. Metody te są dziś w powszechnym użyciu.

Liczby pseudolosowe

W trakcie tego kursu kilkakrotnie wykorzystaliśmy liczby wylosowane przez komputer z zadanego rozkładu. Teraz opiszemy jedną z metod używanych przez współczesne programy komputerowe do ich uzyskiwania. Otrzymane w ten sposób liczby nie są w istocie losowe, lecz są wynikiem działania pewnego dość prostego algorytmu deterministycznego. Dlatego też określa się je często mianem liczb pseudolosowych.

Obecnie wykorzystywane programy komputerowe używają najczęściej następującego algorytmu: dla ustalonych liczb naturalnych a, b i p wybieramy dowolną liczbę naturalną X_0, zwaną ziarnem (ang. seed), a następnie określamy rekurencyjnie ciąg:


X_{n+1} = a X_n + b \pmod{p}.


Mówiąc inaczej, za każdym razem obliczamy X_{n+1}' = a X_n + b, a jako X_{n+1} bierzemy resztę z dzielenia X_{n+1}' przez p - tak więc wszystkie wyrazy naszego ciągu są liczbami naturalnymi mniejszymi od p. Jeżeli parametry a, b, p i X_0 są odpowiednio dobrane, to okazuje się, że liczby:


U_n = X_n /p


mają własności niemal takie same, jak liczby wylosowane z rozkładu jednostajnego na przedziale (0,1).

Istnieją pewne zasady wybierania parametrów. W szczególności, p powinno być bardzo duże, aby jak najbardziej ograniczyć zjawisko okresowości. Z podobnych względów także a powinno być dużą liczbą, najlepiej względnie pierwszą z p. Natomiast wybór b ma mniejsze znaczenie - często przyjmuje się b=0. Okazuje się, że przy odpowiednio wybranych parametrach oraz przy zastosowaniu dodatkowych procedur (patrz poniżej) liczby pseudolosowe i ich zestawy mają bardzo dobre własności - potwierdzają to także odpowiednie testy statystyczne. Przykładowo, program Maple (już w wersji 5) używa generatora liczb pseudolosowych z bardzo dużymi parametrami a oraz p, mianowicie:


a = 427419669081, \;\; p = 999999999989.


Wartość ziarna X_0 można w każdej chwili zadać zgodnie z aktualnymi potrzebami. Może nam na przykład zależeć, aby przy powtórzeniach danego losowania zawsze otrzymywać te same liczby pseudolosowe - zadajemy wtedy taką samą (stałą) wartość X_0 (przy rozpoczęciu nowej sesji z programem Maple wartość ziarna wynosi 1). Z drugiej strony, możemy żądać, aby w każdym losowaniu otrzymywać inne liczby pseudolosowe - można to na przykład osiągnąć, zadając wartość ziarna w zależności od upływu czasu zużytego przez procesor od rozpoczęcia aktualnej sesji.

Bardzo ważne jest też to, aby kolejne liczby pseudolosowe reprezentowały niezależne zmienne losowe. Oczywiście, formalnie one są zależne, jednak przy odpowiednio dobranych parametrach zależność ta jest niezwykle słaba. Co więcej, stosuje się czasem dodatkowe zabezpieczenia polegające na odrzucaniu niektórych liczb.

Mając liczby pseudolosowe z rozkładu jednostajnego na odcinku [0,1], można, stosując odpowiednią metodę, uzyskać liczby pseudolosowe z innego zadanego rozkładu - poniżej opisujemy jedną z takich metod, zwaną metodą odwracania dystrybuanty.

Niech \xi będzie zmienną losową o rozkładzie jednostajnym na przedziale (0,1), czyli:


P(\xi \le x) =\left\{ \begin{array} {ll} 0, &  x < 0\\ x, & 0\le x < 1\\ 1, & 1 \le x, \end{array}  \right.


oraz niech F będzie dystrybuantą interesującego nas w danej chwili rozkładu. Dla uproszczenia załóżmy, że F jest funkcją ściśle rosnącą [AM]. Określamy nową zmienną losową \eta jako:


\eta = F^{-1}(\xi).


Wówczas dla każdego x\in {\Bbb R} mamy:


P(\eta \le x) = P(F^{-1}(\xi) \le x) =  P(\xi  \le  F(x)) = F(x).


Wynika stąd, że F jest dystrybuantą zmiennej losowej \eta. Mówiąc inaczej, jeżeli liczby x_i zostały wylosowane zgodnie z rozkładem jednostajnym na przedziale (0,1), to F^{-1}(x_i) mogą być traktowane jako liczby wylosowane z rozkładu o dystrybuancie F.

Przykład 14.1

Aby wylosować liczbę pseudolosową z rozkładu wykładniczego, rozważamy jego dystrybuantę:


F(x) = 1 - e^{-\lambda x} \ \ \mbox{ dla } x > 0.


W związku z tym, mając daną liczbę pseudolosową u z rozkładu jednostajnego na odcinku (0,1), znajdujemy liczbę pseudolosową:


F^{-1}(u) = - \frac{1}{\lambda}\ln (1- u)


z rozkładu wykładniczego.

Zauważmy, że nie zawsze jest łatwo wyznaczyć efektywnie F^{-1} - jest tak na przykład w przypadku rozkładu normalnego. W takich sytuacjach można szukać przybliżonej wartości F^{-1}(u), rozwiązując numerycznie równanie:


F(x) = u


ze względu na x. Istnieją też inne metody pozyskiwania liczb pseudolosowych z niektórych rozkładów (na przykład rozkładu normalnego), na bazie liczb pochodzących z rozkładu jednostajnego.

Metoda bootstrap

W ostatnich latach stosuje się coraz powszechniej tak zwaną metodę bootstrap. Robi się to wtedy, gdy nie znamy rozkładu, z którego pochodzi próbka, a jej wielkość jest zbyt mała, aby stosować metody oparte na prawach wielkich liczb.

Metoda bootstrap polega na losowaniu kolejno B próbek na podstawie wyjściowej próbki, przy czym losowanie odbywa się ze zwracaniem, a wielkości próbek są takie same jak wielkość próbki wyjściowej. Jeżeli chcemy estymować dany parametr i znamy jego estymator, powiedzmy g(\mbox{X_1, , X_n}), to estymator bootstrapowy danego parametru określamy jako średnią z wartości tego estymatora obliczonych dla każdej próbki, czyli:


\hat{g}(\mbox{x_1, , x_n}) = \frac{1}{B}\sum_{i = 1}^Bg(x^i_1, \dots, x^i_n),


gdzie x^i_1, \dots, x^i_n jest próbką wylosowaną za i-tym razem.

Przykład 14.2

Przypuśćmy, że mamy daną następującą próbkę prostą z nieznanego rozkładu:


2, \ 3, \ 3, \ 3, \ 7, \ 10,\ 12.


Interesuje nas bootstrapowy estymator średniego błędu.

Pamiętamy, że średni błąd wyraża się wzorem:


b = \frac{1}{n}\sum_{i=1}^n|x_i - \bar{x}|.


Losujemy teraz 10 próbek bootstrapowych, przy czym każda taka próbka ma 7 elementów otrzymanych jako wynik losowania ze zwracaniem z podstawowej próbki, jednocześnie wyznaczając dla każdej z nich wartość b. Dostajemy przykładowo:


7,\ 12,\ 12,\ 3,\ 7,\ 10,\ 3; \;\; b\approx 3.10,
3,\ 3,\ 3,\ 10,\ 7,\ 3,\ 2; \;\; b\approx 2.33,
7,\ 3,\ 3,\ 3,\ 3,\ 3,\ 3; \;\; b\approx 0.98,
3,\ 2,\ 3,\ 7,\ 3,\ 7,\ 10; \;\; b\approx 2.57,
3,\ 3,\ 12,\ 3,\ 10,\ 3,\ 2; \;\; b\approx 3.35,
7,\ 3,\ 2,\ 10,\ 7,\ 10,\ 12; \;\; b\approx 2.90,
3,\ 12,\ 3,\ 10,\ 7,\ 12,\ 12; \;\; b\approx 3.51,
3,\ 7,\ 3,\ 3,\ 12,\ 3,\ 3; \;\; b\approx 2.65,
2,\ 2,\ 3,\ 7,\ 7,\ 3,\ 3; \;\; b\approx 1.80,
3,\ 3,\ 3,\ 3,\ 12,\ 3,\ 2; \;\; b\approx 2.24.


Liczymy teraz średnią z tak otrzymanych 10 średnich błędów dla poszczególnych próbek i otrzymujemy estymator bootstrapowy średniego błędu równy w przybliżeniu 2.54.

Na ogół losuje się dużo więcej niż w powyższym przykładzie próbek bootstrapowych - często jest to 1000 próbek.

Metodę bootstrap można używać także do wyznaczania przedziałów ufności określonych parametrów. Istnieje tutaj kilka metod - my poznamy, na przykładzie, tak zwaną metodę percentylową.

Przypuśćmy, że na podstawie naszej próbki chcemy wyznaczyć 90\% przedział ufności dla średniego błędu. Postępujemy wówczas następująco: losujemy 1000 próbek bootstrapowych i dla każdej z nich obliczamy żądany estymator, powiedzmy średni błąd, otrzymując ciąg:


bs_1, \dots, bs_{1000}.


Dla jego zobrazowania rysujemy odpowiedni histogram:



Kwantyle tego ciągu rzędu 0.05 i 0.95 są (z definicji) końcami szukanego przedziału ufności - w naszym przypadku okazał się być nim przedział:
(1.10, 4.08).

Estymacja gęstości

Zdarza się, iż w niektórych sytuacjach nie wiemy z góry, z jakiego rozkładu pochodzi próbka prosta. Warto wówczas oczywiście naszkicować histogram, ale jeżeli próbka wskazuje na to, że rozkład jest ciągły, to warto także sporządzić tak zwany jądrowy estymator gęstości, który w tym punkcie omawiamy. Jego podstawową zaletą (w stosunku do histogramu) jest to, że jest on na ogół funkcją ciągłą.

Funkcję K\colon {\Bbb R} \longrightarrow {\Bbb R} nazywamy jądrem, gdy:

(1) K(x) \ge 0 dla każdego x \in {\Bbb R},

(2) \int_{-\infty}^{\infty} K(x)\,dx = 1,

(3) K(-x) = K(x) dla każdego x \in {\Bbb R}.

Przykładami jąder są następujące funkcje:

  1. K(x) = \frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}} - jądro gaussowskie,
  2. \displaystyle K(x) = \left\{\begin{array} {rl} \frac{3}{4}(1 - x^2), &  |x| < 1\\ 0, & |x| \ge 1 \end{array} \right. - jądro Epanesznikowa,
  3. \displaystyle K(x) = \left\{\begin{array} {rl} 1 - |x|, &  |x| < 1\\ 0, & |x| \ge 1 \end{array} \right. - jądro trójkątne.

Niech będzie dana próbka prosta x_1, \dots, x_nz rozkładu ciągłego o gęstości f. Ustalmy jądro K oraz liczbę h > 0 zwaną dalej szerokością pasma. Funkcję:


\hat{f}(x) = \frac{1}{nh}\sum_{i = 1}^nK\left(\frac{x-x_i}{h}\right)


nazywamy estymatorem jądrowym gęstości f.

Stosując kolejno następujące zmiany zmiennych (patrz wykład z Analizy matematycznej):


t = \frac{x-x_i}{h}\;\; \textrm{dla} i=1,\ldots,n,


można łatwo sprawdzić, że:


\int_{-\infty}^\infty \hat{f}(x)\,dx = 1.


Oczywiście także:


\hat{f}(x)\ge 0\;\;\textrm{dla wszystkich } x \in {\Bbb R},


tak więc funkcja \hat{f} jest gęstością rozkładu prawdopodobieństwa. Zauważmy dodatkowo, że gdy jądro jest funkcją ciągłą, tak jak w podanych powyżej przykładach, to także \hat{f}, jako suma funkcji ciągłych, jest funkcją ciągłą (patrz wykład z Analizy matematycznej).

Można łatwo zrozumieć powód, dla którego jądrowy estymator gęstości \hat{f} powinien dobrze przybliżać daną gęstość f. Otóż, gdy f(x) jest dużą liczbą, to należy się spodziewać, że w pobliżu punktu x znajduje się stosunkowo dużo punktów x_i z naszej próbki, a zatem także wartości K\left(\frac{x-x_i}{h}\right) odpowiadające tym punktom są stosunkowo duże, a więc duża jest też wartość \hat{f}(x). Przeciwnie, tam gdzie wartości f(x) są małe, tam jest niewiele (albo wcale nie ma) punktów próbki i może się zdarzyć, że niemal wszystkie składniki K\left(\frac{x-x_i}{h}\right) będą bliskie zeru, a więc także \hat{f}(x) będzie niewielką liczbą.

Zauważmy jednak, że do wyznaczenia pojedynczej wartości \hat{f}(x) należy, oprócz innych działań, wykonać obliczenia tylu wartości K\left(\frac{x-x_i}{h}\right), ile jest elementów próbki - w praktyce jest to możliwe tylko przy użyciu komputera.


Przykład 14.3

Dla zilustrowania metody estymacji jądrowej wylosujemy 250-elementową próbkę prostą z rozkładu N(50,3), a następnie, wykorzystując jądro gaussowskie, wyznaczymy estymator \hat{f} (kolor czerwony) i narysujemy go na wspólnym rysunku z gęstością rozkładu N(50,3) (kolor niebieski).

Zastosujemy tu trzy różne szerokości pasma h:

1. h = 0.6:



2. h = 1.4:


3. h = 1.05:


Istnieją przesłanki zarówno teoretyczne, jak praktyczne, dotyczące wyboru jądra oraz szerokości pasma, tak aby estymator \hat{f} najlepiej, w pewnym określonym sensie, przybliżał gęstość f. Mówiąc ogólnie, warto jest używać jąder wymienione powyżej, natomiast h można ustalać jako liczbę nieznacznie większą od \frac{\hat{\sigma}}{\sqrt[5]{n}}, gdzie \hat{\sigma} oznacza estymator odchylenia standardowego z próbki, czyli:
\hat{\sigma} = \sqrt{\frac{(x_1-\bar{x}_n)^2 + \dots + (x_n - \bar{x}_n)^2}{n}}.


Warto też zauważyć, że małe wartości h powodują, że wykres \hat{f} jest stosunkowo nieregularny, ale za to "wychwytuje" poszczególne punkty próbki, natomiast dużym wartościom h odpowiada stosunkowo gładki wykres \hat{f}. W powyższym przykładzie \hat{\sigma} \approx 2.92, a więc:


\frac{\hat{\sigma}}{\sqrt[5]{250}} \approx 0.92.