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

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(\mathrm{mod}p).$

Mówiąc inaczej, za każdym razem obliczamy $X_{n^{\prime }+1}=a{X}_n+b$, a jako $X_{n+1}$ bierzemy resztę z dzielenia $X_{n^{\prime }+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)=\{\begin{array}{ll}0,& x<0\\ x,& 0\le x<1\\ 1,& 1\le x,\end{array}$

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 \mathbb{ℝ}$ 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}\text{ 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 Parser nie mógł rozpoznać (błąd składni): {\displaystyle 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:

Parser nie mógł rozpoznać (błąd składni): {\displaystyle \hat{g}(\mbox{x_1, , x_n}) = \frac{1}{B}\sum_{i = 1}^Bg(x^i_1, \dots, x^i_n),}

gdzie $x_1^i,\dots ,x_n^i$ 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-\overline{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\mathrm{\%}$ 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:

$b{s}_1,\dots ,b{s}_{1000}.$

Dla jego zobrazowania rysujemy odpowiedni histogram:

147.eps (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).$