METODY NUMERYCZNE I METODY PROBABILISTYCZNE

12.01.2013

WEKTORY I WŁASNOŚCI WŁASNE

Obliczanie wartości własnych i wektorów własnych macierzy.

Pojęcia podstawowe

RACHUNEK PRAWDOPODOBIEŃSTWA

Test : Prawdopodobieństwo warunkowe i niezależność

WIELKIE UKŁADY RÓWNAŃ LINIOWYCH

Wraz z coraz większymi modelami pojawiającymi się w praktyce obliczeniowej, coraz częściej zachodzi potrzeba rozwiązywania zadań algebry liniowej, w której macierze są co prawda wielkiego wymiaru, ale najczęściej rozrzedzone, to znaczy jest w nich bardzo dużo zer. Bardzo często zdarza się, że macierz wymiaru $\displaystyle N$ ma tylko $\displaystyle O(N)$ niezerowych elementów. Wykorzytanie tej specyficznej własności macierzy nie tylko prowadzi do algorytmów istotnie szybszych od ich analogów dla macierzy gęstych (to znaczy takich, które (w założeniu) mają $\displaystyle N^2$ elementów), ale wręcz są jedynym sposobem na to, by niektóre zadania w ogóle stały się rozwiązywalne przy obecnym stanie techniki obliczeniowej!
Jednym ze szczególnie ważnych źródeł układów równań z macierzami rozrzedzonymi są np. równania różniczkowe cząstkowe (a więc np. modele pogody, naprężeń w konstrukcji samochodu, przenikania kosmetyków do głębszych warstw skóry, itp.).
Modele wielostanowych systemów kolejkowych (np. routera obsługującego wiele komputerów) także prowadzą do gigantycznych układów równań z macierzami rozrzedzonymi o specyficznej strukturze.
Z reguły zadania liniowe wielkiego wymiaru będą miały strukturę macierzy rozrzedzonej, gdyż najczęściej związki pomiędzy niewiadomymi w równaniu nie dotyczą wszystkich, tylko wybranej grupy.

Przykład: Macierz z kolekcji Boeinga

Spójrzmy na macierz sztywności dla modelu silnika lotniczego, wygenerowaną swego czasu w zakładach Boeinga i pochodzącą z dyskretyzacji pewnego równania różniczkowego cząstkowego. Pochodzi z kolekcji Tima Davisa. Jest to mała macierz, wymiaru 8032 (w kolekcji spotkasz równania z milionem i więcej niewiadomych).

Struktura niezerowych elementów macierzy bcsstk38.

Jej współczynnik wypełnienia (to znaczy, stosunek liczby niezerowych do wszystkich elementów macierzy) wynosi jedynie $\displaystyle 0.006$ , a więc macierz jest bardzo rozrzedzona: w każdym wierszu są średnio tylko 44 niezerowe elementy.

PODSTAWY TEORII WERYFIKACJI HIPOTEZ

Formułowanie hipotezy statystycznej rozpoczyna się zebraniem informacji na temat populacji i jej możliwego rozkładu. Dzięki temu możliwe jest zbudowanie zbioru hipotez dopuszczalnych Ω, czyli zbioru rozkładów, które mogą charakteryzować badaną populację. Hipoteza statystyczna to każdy podzbiór zbioru hipotez dopuszczalnych

Podział hipotez statystycznych

Hipotezy statystyczne można podzielić na:

parametryczne - hipoteza dotyczy wartości parametru rozkładu
nieparametryczne - hipoteza dotyczy postaci funkcyjnej rozkładu

Według innego kryterium podział przebiega następująco:

proste - hipoteza jednoznacznie określa rozkład danej populacji, czyli odpowiadający jej podzbiór zbioru Ω zawiera jeden element (rozkład)
złożone - hipoteza określa całą grupę rozkładów, zaś odpowiadający jej podzbiór zbioru Ω zawiera więcej niż jeden element
alternatywna - przyjmujemy ją kiedy odrzucamy hipotezę zerową

Hipotezy weryfikujemy za pomocą testów statystycznych.

Test statystyczny: metoda postępowania, która każdej próbce x₁, x₂, ...,x_n przyporządkowuje ustalonym prawdopodobieństwem decyzje odrzucenia lub przyjęcia sprawdzanej hipotezy.

parametryczne testy istotności - służą do weryfikacji hipotez parametrycznych – odrzucić czy też nie hipotezę wyjściową (zerową);

testy zgodności - testy weryfikujące hipotezy dotyczące zgodności pomiędzy rozkładem wartości w próbce i rozkładem teoretycznym.

Przykładem jest test c² Pearsona;

Parametryczne testy istotności

Rozpatrzymy poniżej 3 parametryczne testy istotności dotyczące:

a) wartości oczekiwanej;

b) różnicy wartości oczekiwanych w dwóch próbkach;

c) wariancji i odchylenia standardowego.

Teza rzeczowa – to, co mamy udowodnić metodą statystyczną, np. że wartość średnia obliczona dla próby jest większa od wartości oczekiwanej w populacji generalnej. W tym celu formułujemy hipotezę, którą zamierzamy weryfikować. Nazywamy ją hipotezą zerową i oznaczamy H₀. Może ona brzmieć następująco: wartość oczekiwana jest równa m₀, co zapiszemy H₀: m=m₀. Zwykle testujemy hipotezę zerową wobec hipotezy alternatywnej H_a, np. H_a: m=m₁¹m₀. Wyniki weryfikacji jakiejś hipotezy nie dają nam absolutnej pewności, ale wnioski możemy sformułować z dowolnie dużym prawdopodobieństwem. Tezę rzeczową, którą chcemy udowodnić metodą statystyczną zwykle nie przyjmujemy jako hipotezy zerowej H_0,ale jako hipotezę alternatywną, którą przyjmujemy po ewentualnym odrzuceniu hipotezy zerowej H₀.

Testowanie składa się z następujących etapów:

1.Sformułowanie tezy rzeczowej i ustaleniu hipotez H₀ i H_a;

2.Wyboru właściwej funkcji testowej (statystyki z próby);

3.Przyjęciu stosownego poziomu istotności a

4.Odczytaniu wartości krytycznych w tablicach dystrybuanty właściwego rozkładu i ustaleniu obszaru krytycznego;

5.Odrzuceniu hipotezy zerowej na korzyść hipotezy alternatywnej, gdy funkcja testowa obliczona z próby znajduje się w obszarze krytycznym i nie odrzucenie jej, gdy funkcja testowa jest poza obszarem krytycznym

Przykład

Należy dokonać oceny partii pudelek zapalek liczącej 100 tys. sztuk. dostawca twierdzi, że w pudelku znajdują się przecietnie 54 zapalki. Zweryfikować hipotezę

H 0 (m = m 0 = 54)

. Ponieważ nie znamy rozkładu liczby zapałek w pudełkach w populacji generalnej, a mozemy łatwo pobrać próbę >= 30 możemy wieć w przybliżeniu skorzystać z rozkładu normalnego. Zkaldamy, że przy próbie o wielkości n = 100 odnotowano średnią arytmetyczną

m n = 51,21

natomiast

σ n' = 2,54

. Weryfikujemy przy poziomie istotności

α = 0,02

ponieważ obraliśmy dużą próbę wiec $\sigma_{M_n} \approx {{\sigma_n'} \over { \sqrt{n}}} = 0{,}245$ . Musimy zatem wyznaczyć t dla ktorego $P \left ( {{|M_n-m_0|} \over {\sigma_{M_n}}} \ge t \right ) = 0{,}02$
Definiujemy unormowaną zmienną Y:
$Y={{M_n-m_0} \over {\sigma_{M_n}}}$
podstawiamy do wzoru
$P(|Y| \ge t) = 0{,}02$
Z własności bezwzględnej wartości:
$P[(Y \ge t) \vee (Y \le -t)] = 0{,}02$
$P(Y \ge t)+P(Y \le -t) = 0{,}02$
Ponieważ funkcja gęstości jest dla rozkładu N(0,1) parzysta to zachodzi równość:
$P(Y \ge t) = P(Y \le -t)$
$2P(Y \ge t) = 0{,}02$
$P(Y \ge t) = 0{,}01$
Wiadomo, że P(A) to to samo co 1-P(A') więc:

1 - P (Y < t) = 0,01

P (Y < t) = 0,99

A P(Y<x) to dystrybuanta - czyli F(x):

F (t) = 0,99

Teraz w tablicy rozkładu normalnego znajdujemy najmniejszą wartość t dla której F(t) wynosi conajmniej 0,99. Jest to wartość 2,33.
Hipotezę

H 0

należy wiec odrzucić na poziomie istotności

α

jeżeli $|M_n - 54| \ge t \cdot 0{,}245$ , w przeciwnym wypadku przy zadanej istotności

α = 0,02

nie możemy ani potwierdzić hipotezy, ani jej odrzucić.
Zgodnie z naszymi danymi wychodzi:

| M n - 54 | = | 51,21 - 54 | = 2,79

$\sigma_{M_n} \cdot t = 0{,}245 \cdot 2{,}33 = 0{,}57085$
Więc:
$| M_n - 54 | \ge \sigma_{M_n} \cdot t$
$2{,}79 \ge 0{,}57085$
Zatem hipotezę możemy odrzucić (jeśli wyjdzie odwrotnie to piszemy że nie odrzucamy ani nie potwierdzamy - tak właśnie trzeba było zrobić na egzaminie, bo wychodziło <).

16.12.2012

KOMPUTEROWE METODY STATYSTYKI

Omówimy 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.

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:

$K(x) = \frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}}$ - jądro gaussowskie,
$\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,
$\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_n$ z 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

$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łą.
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.