Przykład 2.1 (Rozkład prawdopodobieństwa i próbka losowa)

Wylosujmy ,,próbkę” średniego rozmiaru, powiedzmy n=100 z rozkładu normalnego N⁢0,1:

> n <- 100
> X <- rnorm(n)
> X

(aby przyjrzeć się funkcji rnorm, napiszmy ?rnorm). Wektor X zawiera ,,realizacje” niezależnych zmiennych losowych X1,…,X100. Możemy obliczyć średnią X¯=n-1⁢∑i=1nXi i wariancję próbkową X¯=n-1-1⁢∑i=1nXi-X¯2. Jakich wyników oczekujemy? Zobaczmy:

> mean(X)
> var(X)

Porównajmy kwantyl empiryczny rzędu, powiedzmy p=1/3, z teoretycznym:

> quantile(X,p)
> qnorm(p)

Możemy powtórzyć nasze ,,doświadczenie losowe” i zobaczyć jakie są losowe fluktuacje wyników. Warto wykonać nasz mini-programik jeszcze raz, albo parę razy. Teraz sprobujmy ,,zobaczyć” próbkę. Najlepiej tak:

> hist(X,prob=TRUE) ( proszę się dowiedzieć co znaczy parametr `prob' ! )
> rug(X)

Możemy łatwo narysować wykres gęstości. Funkcja obliczająca gęstość φ rozkładu N⁢0,1 nazywa się dnorm, a curve jest funkcją rysującą wykresy.

> curve(dnorm(x),col="blue",add=TRUE)

(nawiasem mówiąc, zachęcam początkujących pRobabilistów do oznaczania wygenerowanych wektorów losowych dużymi literami, np. X, a deterministycznych zmiennych i wektorów – małymi, np. x, podobnie jak na rachunku prawdopodobieństwa).

Podobnie, możemy porównać dystrybuantę empiryczną z prawdziwą dystrybuantą Φ. Funkcja obliczająca dystrybuantę empiryczną nazywa się ecdf, zaś dystrybuantę rozkładu N⁢0,1 – pnorm.

> plot(ecdf(X))
> curve(pnorm(x),from=xmin,to=xmax,col="blue",add=TRUE)

Test Kołmogorowa-Smirnowa oblicza maksymalną odległość D=∑xF⌃n⁢x-F⁢x, gdzie F⌃n jest dystrybuantą empiryczną i F=Φ.

> ks.test(X,pnorm,exact=TRUE)

Przypomnijmy sobie z wykładu ze statystyki, co znaczy podana przez test p-wartość. Jakiego wyniku ,,spodziewaliśmy się”?

Przykład 2.2 (Powtarzanie doświadczenia symulacyjnego)

Symulacje polegają na powtórzeniu całego doświadczenia wiele razy, powiedzmy m=10000 razy, zanotowaniu wyników i uznaniu powstałego rozkładu empirycznego za przybliżenie badanego rozkładu prawdopodobieństwa. Bardzo ważne jest zrozumienie różnej roli jaką pełni tu n (rozmiar próbki w pojedynczym doświadczeniu, a więc ,,parametr badanego zjawiska”) i m (liczba powtórzeń, która powinna być możliwie największa aby zwiększyć dokładność badania symulacyjnego). Następujący progarmik podkreśla logiczną konstrukcję powtarzanego doświadczenia:

> m <- 10000
> n <- 100
>
> # Przygotowujemy wektory w którym zapiszemy wyniki:
> D <- c()
> P <- c()
>
> for (i in 1:m)
> {
> X <- rnorm(n)
> Test <- ks.test(X,pnorm,exact=TRUE)
> D[i] <- Test$statistic
> P[i] <- Test$p.value
> } # koniec pętli for
>
> # Analizujemy wyniki:
> hist(D, prob=TRUE)
> hist(P, prob=TRUE)

Co prawda powyższy programik spełnia swoją rolę, jednak struktura pakiektu R jest dostosowana do innego stylu pisania programów. Zamiast pętli for zalecane jest używanie funkcji które powtarzają pewne operacje i posługiwanie się, kiedykolwiek to możliwe, całymi wektorami a nie pojedynczymi komponentami. Poniższy fragment kodu zastępuje pętlę for

> DiP <- replicate(m, ks.test(rnorm(n),pnorm,exact=TRUE)[1:2]))
>
> DiP <- t(DiP) # transpozycja macierzy ułatwia oglądanie wyników
> D <- as.numeric((DiP)[,1]) # pierwsza kolumna macierzy
> P <- as.numeric((DiP)[,2]) # druga kolumna macierzy
> # obiekty typu ,,list'' przerobimy na wektory liczbowe:
> D <- as.numeric(D); P <- as.numeric(D)

Przykład 2.3 (Mocne Prawo Wielkich Liczb)

Wylosujmy próbkę z ,,jakiegoś” rozkładu prawdopodobieństwa. Weźmy na przykład niezależne zmienne o rozkładzie wykładniczym, X1,…,Xn,…⁢Ex⁢2. Obliczmy średnie empiryczne

Mn=Snn=1n⁢∑i=1nXi.

Zróbmy wykres ciągu średnich S1/1,S2/2,…,Sn/n,….

> nmax <- 1000 # komputerowy odpowiednik ,,n→∞''
> n <- (1:nmax)
> lambda <- 2 > X <- rexp(nmax,rate=lambda)
> S <- cumsum(X) # ciąg narastających sum
> M <- S/n # działania w R (np. dzielenie) są wykonywane ,,po współrzędnych''
> plot(n,M,type="l") # ,,zaklęcie'' type="l" powoduje narysowanie łamanej

Teraz spóbujmy podobne doświadczenie zrobić dla zmiennych Xi=1/Ui-1, gdzie Ui∼U⁢0,1 (Xi są próbką z tak zwanego rozkładu Pareto).

> # potrzebujemy dużej próbki, żeby się zorientować, co się dzieje...
> nmax <- 100000
> n <- (1:nmax)
> X <- 1/runif(nmax)-1
> M <- cumsum(X)/n
> plot(log(n),M,type="l") # i zrobimy wykres w skali logarytmicznej

Przykład 2.4 (Centralne Twierdzenie Graniczne)

CTG jest jednym ze ,,słabych” twierdzeń granicznych rachunku prawdopodobieństwa, to znaczy dotyczy zbieżności rozkładów. Symulacyne ,,sprawdzenie” lub ilustracja takich twierdzeń wymaga powtarzania doświadczenia wiele razy, podobnie jak w Przykładzie 2.2. Pojedyncze doświadczenie polega na obliczeniu sumy

Sn=∑i=1nXi,

gdzie X1,…,Xn jest próbką z ,,jakiegoś” rozkładu prawdopodobieństwa i n jest ,,duże”. Weźmy na przykład niezależne zmienne o rozkładzie wykładniczym, jak w Przykładzie 2.3.

> m <- 10000
> n <- 100
> lambda <- 2
> S <- replicate(m, sum(rexp(n,rate=lambda)))
> hist(S,prob=TRUE)
> curve(dnorm(x,mean=n/lambda,sd=sqrt(n)/lambda),col="blue",add=TRUE)
> # wydaje się na podstawie obrazka, że dopasowanie jest znakomite
> ks.test(S,pnorm,mean=n/lambda,sd=sqrt(n)/lambda)
> # ale test Kołmogorowa-Smirnowa ,,widzi'' pewne odchylenie od rozkładu normalnego >