1 Generatory liczb losowych w R

W wielu zagadnieniach związanych z modelowaniem metodami Monte Carlo złożonych zjawisk (np. finansowych, biologicznych, fizycznych, chemicznych, metrologicznych, itp.) potrzebne są narzędzia umożliwiające symulować losowość, tak aby odzwierciedlić stochastyczną naturę zjawiska.

Generowane przez komputer liczby nazywamy pseudolosowymi, ponieważ mają one symulować losowość, jednakże są one wyznaczane w sposób deterministyczny (choć często w dość skomplikowany sposób). Ponieważ zapotrzebowanie na na dobre liczby pseudolosowe nieustannie rośnie w literaturze można znaleźć wiele metod ich generowania.

1.1 Liniowa metoda kongruencji

Mimo, że będziemy głównie korzystać z metod zaimplementowanych w R, przedstawimy poniżej jeden z najbardziej klasycznych algorytmów do symulacji liczb losowych zwany liniową metodą kongruencji (linear congruential method). Główna zaletą liniowych generatorów kongruencjalnych jest ich stosunkowo łatwą implementację i szybkość oraz zapotrzebowanie na niewielką ilość pamięci. Warto jednak zaznaczyć, że inne metody, takie jak Mersenne Twister, są dziś znacznie bardziej powszechne w praktycznym użyciu.

Liniowe generatory kongruencyjne są definiowane przez relację rekurencyjną: \[X_{i+1}=(aX_{i}+c)\text{ mod }m,\] gdzie \(a, c, m \in \mathbb{N}_{+}\) są parametrami, które należy wybrać oraz wartość \(X_0\) będącą tak zwanym ziarnem (seed) algorytmu.

Zatem, algorytm generuje liczby z przedziału \([0,m-1]\).

1.2 Zadanie 1 - kongruencyjny generator liczb losowych

Zaimplementuj funkcję lcm(N, x0, a, c, m) obliczającą według powyższego wzoru liniowy generator kongurencyjny, gdzie \(N\) jest długością generowanego ciągu, a pozostałe argumenty parametrami generatora. Przeskaluj wartość \(X_n\) tak, aby funkcja zwracała wartość z przedziału \([0,1]\). Wywołaj zaimplementowaną funkcję dla podanych wartości argumentów lcm(N = 10, x0 = 4, a = 13, c = 0, m = 64), a następnie narysuj otrzymane wyniki. Co można zauważyć na podstawie wykresu? Czy wygenerowane liczby są według Ciebie losowe?

Uwaga: Teoria dotycząca optymalnego wyboru parametrów i liczby początkowej \(X_0\) jest dość skomplikowana. Powszechnym wyborem jest jednak jest przyjęcie za \(X_0\) bieżącego czasu systemowego w mikrosekundach.

1.3 Zadanie 2 - kongruencyjny generator liczb losowych

Wywołaj wcześniej napisana funkcję dla poniższych argumentów. Czy w tym przypadku otrzymane wartości są bardziej losowe?

  lcm(N = 10, x0 = as.numeric(Sys.time()), a = 1103515245, c = 12345, m = 2 ** 32)

Uwaga: W ten sposób możemy generować losowe wartości z przedziału jednostajnego \(U[0,1]\).

2 Popularne rozkłady zmiennych.

W języku R jest wiele funkcji do obsługi większości rozkładów zmiennych losowych. Nazewnictwo funkcji jest zestandaryzowane i można je opisać za pomocą poniższego schematu:

[prefix][nazwa_rodziny_rozkładów]()

Możliwe wartości dla prefixów to:

Pozostałe argumenty określają parametry rozkładu w wybranej rodzinie rozkładów.

2.1 Rozkład jednostajny

Zacznijmy od rozkładu jednostajnego, który jest jednym z najprostszych i najbardziej podstawowych rozkładów. Posiada on poniższą funkcję gęstości: \[f(x)= \begin{cases}\frac{1}{b-a}, & \text { gdy } a \leq x \leq b \\ 0, & \text { gdy } x<a \text { lub } x>b.\end{cases}\]

unif jest nazwą rodziny rozkładów dla rozkładów jednostajnych.

Na przykład aby wygenerować 10 zmiennych losowych z rozkładu \(U[0,1]\) wystarczy:

set.seed(5)
runif(10) # wartości parametrów w tym przypadku są domyślne, tj. a=0, b=1
##  [1] 0.2002145 0.6852186 0.9168758 0.2843995 0.1046501 0.7010575 0.5279600
##  [8] 0.8079352 0.9565001 0.1104530

Funkcja set.seed ustawia ziarno i determinuje wartości kolejnych liczb pseudo losowych. Dzięki temu możemy zapewnić możliwość powtórzenia swoich wyników, które bazują na losowości.

Dla rozkładu jednostajnego na innym przedziale niż \([0{,}1]\), używamy:

runif(10, min=2, max=8) # 1 wersja
runif(10,2,8)           # 2 wersja

2.2 Zadania - rozkład jednostajny

2.2.1 Zadanie unif - 4 wykresy: gęstości i dystrybuanty

Stwórz cztery wykresy na wspólnym rysunku (używając par(mfrow=c(2,2)) do podziału 1 rysunku na 4 panele) ukazujące dwie gęstość rozkładu jednostajnego dla różnych parametrów i odpowiednio ich dystrybuanty.

Wynik dla wybranych parametrów powinien przypominać

2.2.2 Zadanie unif - funkcja kwantylowa

Narysuj funkcję kwantylową dla zmiennych z przedziału [0,1].

Funkcja kwantylowa jest odwrotnością dystrybuanty. Dla danego prawdopodobieństwa \(p\) wyznacza wartość \(q\) dla której prawdopodobieństwo, że zmienna losowa osiągnie wartość mniejszą lub równą \(q\) wynosi \(p\), to jest, \(P(X\leq q)=p\).

2.2.3 Zadanie unif - histogram i gęstość rozkładu

Porównaj histogram z prawdziwą gęstością rozkładu. Wygeneruj 100 obserwacji z rozkładu \(U[2,8]\) i przedstaw je jako histogram. Następnie oblicz i nanieś na niego teoretyczny wykres gęstości rozkładu \(U[2,8]\).

Rysunek powinien przypominać (czerwona linia to teoretyczna gęstość uzyskana za pomocą abline. Histogram, by zamiast liczności umieścił na wykresie częstotliwość występowania, powienien mieć argument probability = TRUE):

2.3 Metoda dystrybuanty odwrotnej

Najprostszą metodą symulacji obserwacji z ogólnych zmiennych losowych jest tak zwana metoda dystrybuanty odwrotnej.

Załóżmy, że interesuje nas zmienna losowa \(X\), której posiada dystrybuantę \(F_X\). Musimy założyć, że dystrybuanta \(F_X\) jest:

  1. znana,

  2. ciągła,

  3. ściśle rosnąca,

wtedy \[X=F^{-1}_X(U),\]

gdzie \(F_X^{-1}\) jest odwrotnością \(F_X\) oraz \(U\) jest zmienną losową z rozkładu jednostajnego (z ang. uniform) \(U[0,1]\) tzn. (\(U\sim U[0,1]\)).

Intuicja

Szukamy ściśle monotonicznego odwzorowania \(T\colon [0,1]\rightarrow\mathbb{R}\), takiego, że zmienna \(X\) ma rozkład \(T(U)\), gdzie \(U\) ma rozkład jednostajny \(U[0,1]\). Zatem, \[F_X(x)=P(X\leq x)=P\left(T(U)\leq x\right)=P\left(U\leq T^{-1}(x)\right)=T^{-1}(x),\] gdzie ostatnia równość wynika z faktu, że \(P(U\leq y)=y\), gdy \(U\sim U[0,1]\).

Stąd wynika, że \(T(u)=F^{-1}_X(u)\).

W rezultacie mamy następujący algorytm na generowanie zmiennych losowych \(X\) o rozkładzie \(F_X\):

  1. oblicz \(F^{-1}_X\) korzystając z \(F\),

  2. wygeneruj niezależny zmienne losowe \(u_1, \dots, u_n\) z rozkładu jednostajnego \(U[0,1]\),

  3. oblicz \(x_1=F^{-1}_X(u_1), \dots, x_n=F^{-1}_X(u_n)\).

Wtedy zmienne losowe \(x_1, \dots, x_n\) są niezależne o rozkładzie \(F_X\).

2.3.1 Zadanie - metoda dystrybuanty odwrotnej + rozkład wykładniczy

Wygeneruj 5 obserwacji z rozkładu wykładniczego metodą dystrybuanty odwrotnej. Przypomnijmy, że jeśli \(X\) ma rozkład wykładniczy z parametrem \(\lambda\), to dystrybuanta wyraża się poniższym wzorem \[F_X(x)=1-e^{-\lambda x}, \quad \text{dla }x\geq 0.\]

Wskazówka: należy samodzielnie wyznaczyć funkcję \(F^{-1}_{X; \lambda}(u)\).

Przykładowy wynik (przy ustalonym set.seed(202) i dla lambda=2) to

## [1] 0.06897777 0.55456051 0.20008199 0.11390671 0.22965812

2.3.2 Metoda dystrybuanty odwrotnej cd

Wiemy, jak generować rozkład jednostajny na przedziale \([0,1]\). Pytanie jak można teraz zadać, to jak wygenerować przedział jednostajny na przedziale ogólnym \([a,b]\).

Możemy również użyć metody dystrybuanty odwrotnej. Przypomnijmy, że dystrybuanta rozkładu jednostajnego wyraża się wzorem: \[F(x)=\frac{x-a}{b-a}.\] Zatem, \[F^{-1}(u)=a+(b-a)u.\] Wobec tego mając zmienne \(u_1, \dots, u_n\) z rozkładu jednostajnego na przedziale \([0,1]\) możemy symulować rozkład \(U[a,b]\) za pomocą \[x_i=a+(b-a)u_i\]

2.3.3 Zadanie - metoda dystrybuanty odwrotnej - rozkład jednostajny

Używając metody dystrybuanty odwrotnej wygeneruj \(5\) zmiennych z rozkładu jednostajnego na przedziale \([2,6]\). Następnie wygeneruj \(5\) obserwacji z rozkładu \(U[2,6]\) używając gotowej funkcji w R.

Przykładowy wynik przy set.seed(202):

## [1] 2.515450 4.680605 3.319159 2.814909 3.473138
## [1] 5.129214 2.793936 2.876312 3.767799 4.569107

2.3.4 Uwagi - metoda dystrybuanty odwrotnej

Większość generatorów liczb losowych na komputerze, generuje tylko liczby z rozkładu jednostajnego, a następnie korzystając z metody dystrybuanty odwrotnej oblicza zmienne z innych rozkładów.

Wyjątkiem jest rozkład normalny, który nie ma prostej postaci analitycznej rozkładu, więc aby użyć metody dystrybuanty odwrotnej trzeba wykorzystać wymagające obliczeniowo aproksymacje numeryczne.

2.4 Rozkład normalny

Jednym z ważniejszych rozkładów jest rozkład normalny (Gaussowski) \(N(\mu, \sigma^2)\). Ma on dwa parametry: średnią \(\mu\) oraz wariancję \(\sigma^2\).

Uwaga: W języku R, funkcje dla rozkładu normalnego przyjmują jako parametr rozkładu odchylenie standardowe, czyli pierwiastek kwadratowy z wariancji, a nie samą wariancję.

Jeżeli \(X\sim N(0,1)\), to mówimy, że zmienna losowa \(X\) ma rozkład standardowy normalny. Funkcja gęstości rozkładu normalnego ma następujący wzór: \[f(x)=\frac{1}{\sigma \sqrt{2\pi}}e^{-\frac{(x-\mu)^2}{2\sigma^2}}.\]

Właściwości:

  • Symetria - rozkład normalny jest symetryczny względem średniej \(\mu\).

  • Krzywa rozkładu normalnego ma charakterystyczny dzwonowy kształt.

  • Reguła 68-95-99.7 - około 68% wartości mieści się w przedziale \([\mu-\sigma, \mu+\sigma]\), 95% w przedziale \([\mu-2\sigma, \mu+2\sigma]\), a 99.7% w przedziale \([\mu-3\sigma, \mu+3\sigma]\).

Zastosowania:

  • testy statystyczne, np. ANOVA, t-Studenta,

  • przedziału ufności,

  • analiza regresji,

  • modele predykcyjne.

Sprawdźmy, jak wyglądają podstawowe funkcje dotyczące rozkładu normalnego.

par(mfrow=c(2,2))
mu=0.5
sigma=1
x <- seq(-4,4,0.1)
values <- seq(0,1,0.1)
plot(x, dnorm(x, mean = mu, sd = sigma), type="l",  main = "Gęstość") # domyślnie oblicza dla rozkładu standardowego normalnego 
plot(x,pnorm(x, mean = mu, sd = sigma), type="l", main = "Dystrybuanta")
plot(values, qnorm(values, mean = mu, sd = sigma), type="l", main = "Funkcja kwantylowa")

draws <- rnorm(100, mean = mu, sd = sigma)
hist(draws, probability = TRUE)
curve(dnorm(x, mean = mu, sd = sigma), add = TRUE, col = "red")

par(mfrow=c(1,1))

2.4.1 Zadania - rozkład normaly

2.4.1.1 Zadanie - wykresy rozkładu normalnego

Przetestuj powyższy kod dla innych parametrów rozkładu normalnego. Jak zmiana parametrów wpływa na kształt wykresów?

2.4.1.2 Zadanie - symulacja egzaminu

Załóżmy, że wyniki egzaminu przeprowadzonego dla \(200\) studentów, mają rozkład normalny o średniej \(72\) i odchyleniu standardowym \(15.2\).

  1. Wygeneruj wyniki z egzaminu.

  2. Wyznacz, jakie jest prawdopodobieństwo uzyskania \(84\) punkty lub więcej na egzaminie. Wskazówka z rachunku prawdopodobieństwa $P(x ) = $

## [1] 0.2149176

2.4.1.3 Zadanie - 3 odchylenia standardowe

Dla zmiennej losowej mającej rozkład normalny \(N(0.5, 4)\) wyznacz prawdopodobieństwo, że znajduje się ona w granicach \(3\) odchyleń standardowych od średniej (tj. wyznacz \(P(\mu-3\sigma<X< \mu+3\sigma)\)).

Wskazówka:

\(P(\mu-3\sigma<X< \mu+3\sigma)=\dots\)

## [1] 0.9973002

2.5 Rozkład chi-kwadrat \(\chi^2\)

Gdy \(X_1, \dots, X_m\) są niezależnymi zmiennymi losowymi o rozkładzie \(N(0,1)\), wtedy \[V=X_1^2+\dots + X_m^2\] ma rozkład chi-kwadrat o \(m\) stopniach swobody \(\chi^2_{(m)}\).

Właściwości:

  • Rozkład chi-kwadrat jest asymetryczny i ma jedno ramię.

  • Jego kształt zależy od liczby stopni swobody. Przy małej liczbie stopni swobody jest bardziej skośny, a przy większej liczbie stopni swobody staje się bardziej zbliżony do rozkładu normalnego.

Zastosowania:

  • Testy zgodności, służą do sprawdzania, czy obserwowane dane pasują do oczekiwanego rozkładu teoretycznego.

  • Testy niezależności, służą do badania, czy istnieje zależność między dwiema zmiennymi kategorycznymi.

  • ANOVA czyli analiza wariancji w porównaniu istotności grup jakościowych.

Nazwa rodziny rozkładów chi-kwadrat w R to chisq.

2.5.1 Zadanie - rozkład chi-kwadrat

Narysuj wykres gęstości rozkładu chi-kwadrat o 7 stopniach swobody. Następnie porównaj go z wykresem chi-kwadrat o 2 i 20 stopniach swobody.

Wynik powinien przypominać:

2.6 Rozkład t-Studenta

Jest również nazywany rozkładem t. Bardzo użyteczny przy pracy z małymi próbkami. Jest on podobny kształtem do rozkładu normalnego, natomiast ma tzw. “grubsze ogony”. Oznacza to, że prawdopodobieństwo otrzymania skrajnych obserwacji jest większe niż w przypadku rozkładu normalnego.

Zmienna T \[T=\frac{Z}{\sqrt{V/m}}\] ma rozkład t-Studenta o \(m\) stopniach swobody \(t_{(m)}\), gdzie \(Z\) jest zmienną losowa o rozkładzie \(N(0,1)\), \(V\) ma rozkład \(\chi^2_{(m)}\) oraz \(Z\) i \(V\) są niezależne.

W języku R rodzinę rozkładów t-Studenta wywołuje się za pomocą t.

Uwaga! Jeżeli próba prosta \(X_1, \dots, X_n\) ma rozkład \(N(\mu, \sigma^2)\), wtedy \[t=\frac{\bar{x}-\mu}{\sqrt{s^2/n}}\sim t_{(n-1)},\] gdzie \(s^2\) jest estymatorem wariancji \(s^2=\frac{1}{n-1}\sum_{i=1}^n \left(X_i-\bar{X}\right)^2\).

Zatem rozkład t-Studenta w praktyce często się wykorzystuje, gdy wariancja jest nieznana i musimy użyć estymatora i za jego pomocą zestandaryzować dane.

Właściwości:

  • średnia dla rozkładu t-Studenta wynosi \(0\),

  • liczba stopni swobody \(m\) determinuje kształt rozkładu. Im większa liczba stopni swobody, tym bardziej rozkład t-Studenta przypomina rozkład normalny.

Zastosowania:

  • Test t-Studenta dla jednej próby - porównuje średnią próby z wartością populacji, gdy liczba obserwacji jest mała, a populacja ma rozkład normalny lub zbliżony do normalnego.

  • Używany do porównywania średnich z dwóch niezależnych prób.

  • Używany do budowy przedziałów ufności, gdy liczba obserwacji jest mała, a odchylenie standardowe populacji nie jest znane.

2.6.1 Zadania - rozkład t-Studenta

2.6.1.1 Zadanie - wykresy rozkładu t-Studenta

Narysuj wykres rozkładu standardowego normalnego, a następnie nanieś na niego rozkład t-Studenta o \(2\), \(5\) i \(10\) stopniach swobody. Porównaj otrzymane wykresy.

2.6.1.2 Zadanie - statystyka t

Wygeneruj \(10\) obserwacji z rozkładu normalnego o średniej \(0.5\) i wariancji \(2\), a następnie oblicz statystykę \(t\) (ze wzoru powyżej). Narysuj rozkład t-Studenta o \(9\) stopniach swobody i zaznacz obliczoną wartość statystyki \(t\) na wykresie.

Na poniższym wykresie statystyka \(t\) została zaznaczona pionową, czerwoną linią (polecenie abline):

2.7 Rozkład F

Nazywany jest również rozkładem F - Snedecora lub rozkład Fishera-Snedecora.

Jeżeli \(V_1\) i \(V_2\) są dwiema niezależnymi zmiennymi losowymi o rozkładzie \(\chi^2\) o stopniach swobody \(m_1\) oraz \(m_2\) odpowiednio, wtedy \[F=\frac{V_1/m_1}{V_2/m_2}\] ma F rozkład \(F_{(m_1, m_2)}\).

Właściwości:

  • rozkład F jest niesymetryczny i ma ogon skierowany w prawo,

  • nieujemny.

Zastosowania:

  • ANOVA,

  • testowanie hipotez dotyczących proporcji wariancji dwóch populacji.

Funkcja gęstości i dystrybuanty dla rozkładu F w zależności od różnych wartosci parametrów prezentuje się następująco.

x=seq(0, 8, 0.01)
par(mfrow=c(1,2))
plot(x, df(x, 5,3), type="l", main ="Gęstość F(5,3)")
lines(x, df(x, 5,1), type="l", main ="Gęstość F(5,1)", col="red")
plot(x, pf(x, 5,3), type="l", main ="Dystrybuanta F(5,3)")
lines(x, pf(x, 5,1), type="l", main ="Dystrybuanta F(5,1)", col="red")

par(mfrow=c(1,1))

2.8 Rozkład dwumianowy (rozkład Bernoulliego)

Jest to rozkład dyskretny. Opisuje wynik \(n\) niezależnych prób w eksperymencie. Zakłada się, że każda próba ma tylko dwa wyniki, tj. sukces lub porażkę. Jeśli prawdopodobieństwo sukcesu wynosi \(p\), to prawdopodobieństwo uzyskania \(k\) udanych wyników w eksperymencie składającym się z \(n\) niezależnych prób jest następujące: \[P(X=k)=\binom{n}{k}p^k(1-p)^{n-k}.\] Rodzina rozkładów Bernoullego w języku nazywa się binom.

Wygenerowanie zmiennych z rozkładu Bernoulliego przy pomocy rozkładu jednostajnego \(U[0,1]\) jest również możliwe.

Algorytm:

  1. Generujemy liczby \(u_1,\dots, u_n\) z rozkładu \(U[0,1]\).

  2. Następnie przekształcamy wylosowane liczby na zmienne binarne. Jeżeli \(u_i<p\) to przypisujemy wartość \(0\), a w przeciwnym przypadku, przypisujemy wartość \(1\).

2.8.1 Zadania - rozkład Bernoullego

2.8.1.1 Zadanie - quiz z angielskiego

W quizie z języka angielskiego jest dwanaście pytań wielokrotnego wyboru. Każde pytanie ma pięć możliwych odpowiedzi i tylko jedna z nich jest poprawna. Znajdź prawdopodobieństwo udzielenia czterech lub mniej poprawnych odpowiedzi, jeśli uczeń spróbuje odpowiedzieć na każde pytanie losowo.

Wskazówka: Prawdopodobieństwo udzielenia poprawnej odpowiedzi na losowo wybrane pytanie wynosi \(1/5\). Aby znaleźć prawdopodobieństwo uzyskania czterech lub mniej poprawnych odpowiedzi w losowych próbach, stosujemy funkcję dbinom z \(x = 0, \dots,4\).

## [1] 0.9274445

2.8.1.2 Zadanie - Bernulli - frakcje

Utwórz wektor \(x\) o rozmiarze \(50\), zawierający losowe liczby z przedziału \([0,1]\). Następnie utwórz nowy wektor \(y\) o tym samym rozmiarze, w którym zapiszesz \(0\) jeżeli wartość \(x_i<0.5\) i \(1\) w przeciwnym przypadku. Oblicz frakcję jedynek w wektorze \(y\). Następnie wygeneruj wektor \(z\) używając funkcji rbinorm o prawdopodobieństwie sukcesu równemu \(\frac{1}{2}\) i również dla niego oblicz frakcję otrzymanych jedynek. Czy wyniki są podobne? Wykonaj całe zadanie bez używania pętli. Zwiększ rozmiar wektorów do 500 i powtórz symulację.

Przykładowy wynik (dla długości równej 50):

## [1] 0.42
## [1] 0.56

2.9 Rozkład Poissona

Jest to drugi przykład rozkładu dyskretnego. Przypomnijmy, że rozkład Bernoulliego posiada górna granice (\(n\)) to jest liczbę eksperymentów. Czasami przeprowadzamy eksperyment i liczymy sukcesy bez zadanego z góry ograniczenia. Eksperymenty te są zwykle przeprowadzane przez pewien czas lub na pewnej przestrzeni, lub w obu przypadkach, na przykład:

  • liczba zdjęć zaobserwowanych przez detektor w ciągu minuty,

  • liczba włączeń klimatyzacji w ciągu godziny,

  • liczba budynków w mili kwadratowej,

  • liczba tranzystorów na płytce drukowanej.

Rozkład Poissona opisuje rozkład prawdopodobieństwa wystąpienia niezależnych zdarzeń w przedziale. Jeżeli parametr \(\lambda\) jest średnim wystąpieniem na przedział, to prawdopodobieństwo wystąpienia \(k\) zdarzeń w danym przedziale wynosi: \[P(X=k)=\frac{\lambda^k e^{-\lambda}}{k!}.\] Właściwości:

  • Rozkład Poissona jest dyskretny i definiowany na nieujemnych liczbach całkowitych: \(0,1,2,\dots\).

  • Wszystkie momenty wyższego rzędu rozkładu Poissona są równe \(\lambda\).

  • Własność braku pamięci, tj. prawdopodobieństwo wystąpienia zdarzenia w przyszłości jest niezależne od przeszłości.

Zastosowania:

  • modelowania liczby zdarzeń,

  • analizy danych, gdzie liczba zdarzeń jest niewielka w stosunku do wielkości badanej populacji.

Spójrzmy na dystrybuantę rozkładu Poissona.

plot(x, ppois(x, lambda = 2), type="s", main = "Dystrybuanta dla rozkładu Poissona(2)")

2.9.1 Zadanie - Poisson - pożary

W mieście, średnia liczba zgłoszeń pożarów wynosi 3 dziennie. Zakładając, że liczba zgłoszeń pożarów ma rozkład Poissona:

  • Oblicz prawdopodobieństwo, że w ciągu jednego dnia zostanie zgłoszonych dokładnie 5 pożarów.

  • Jaka jest szansa, że nie zostanie zgłoszony żaden pożar w ciągu jednego dnia?

  • Oblicz prawdopodobieństwo, że w ciągu tygodnia (7 dni) liczba zgłoszeń pożarów będzie większa niż 25.

Wskazówka: Dla rozkładu Poissona zachodzi tak zwana zasada addytywności, to znaczy jeżeli mamy dwa rozkłady Poisson(\(\lambda_1\)) i Poisson(\(\lambda_2\)), to ich suma ma rozkład Poissona z parametrem \(\lambda_1+\lambda_2\).

  • Wygeneruj 1000 symulacji liczby zgłoszeń pożarów na tydzień (7 dni). Narysuj histogram wyników i porównaj go z teoretycznym rozkładem Poissona o parametrach odpowiadających tygodniowi.

Przykładowy histogram:

3 Analiza rozkładu danych

Znalezienie odpowiedniego rozkładu dla danych jest kluczowe w statystyce, ponieważ umożliwia poprawne modelowanie, analizę i przewidywanie. Na przykład można go użyć do:

Dla danego zestawu danych powinno się wykonać następujące kroki:

  1. Znalezienie kandydata na rozkład.

  2. Estymacja parametrów rozkładu.

  3. Ocena dopasowania rozkładu.

Powyższy schemat opiszemy bardziej szczegółowo w poniższych podsekcjach.

Dopasowany rozkład danych możemy wykorzystanie na przykład do analizy lub tworzenia predykcji.

3.1 Znalezienie kandydata na rozkład

Istotna jest eksploracja danych i ich rozumienie. W tym celu możemy wykorzystać statystyki opisowe i wizualizacje (opisane w laboratorium 6), jak na przykład funkcja summary(), histogramy lub wykresy pudełkowe.

Funkcja descdist() z pakietu fitdistrplus jest bardzo użytecznym narzędziem do dopasowania odpowiedniego rozkładu do danych. Funkcja ta generuje wykres i wypisuje informacje, które pomagają w ocenie, jak dobrze dane empiryczne pasują do wybranych rozkładów teoretycznych.

3.2 Estymacja parametrów rozkładów

W praktyce często nie generujemy sami rozkładu ze znanymi parametrami, a raczej otrzymujemy obserwacje i chcemy dopasować najlepszy rozkład z najlepiej dobranymi parametrami. Wyznaczanie nieznanych parametrów rozkładów statystycznych odbywa się za pomocą różnych metod estymacji. Przedstawimy kilka najczęściej stosowanych metod.

3.2.1 Metoda największej wiarygodności (MLE)

Polega ona na znalezieniu wartości parametrów, które maksymalizują funkcję wiarygodności, tj. prawdopodobieństwo uzyskania obserwowanych danych przy danych parametrach.

W R możemy użyć funkcji fitdist() z biblioteki fitdistrplus. Pierwszy argumentem funkcji fitdistr() jest wektor obserwacji, na bazie którego wykonana będzie estymacja, a drugim powinna być nazwa rozkładu. Trzecim argumentem jest nazwa metody jaką chcemy użyć do estymacji parametru. Domyślnie jest to metoda największej wiarygodności. W rezultacie otrzymamy obiekt klasy fitdist, który zawiera estymowane parametry rozkładu (estimate) i błąd standardowy dla tych parametrów (Std. Error).

library(fitdistrplus)
wek <- rnorm(100,0.2,3)
fitdist(wek, "norm")
## Fitting of the distribution ' norm ' by maximum likelihood 
## Parameters:
##       estimate Std. Error
## mean 0.4419066  0.2857308
## sd   2.8573081  0.2020421

3.2.2 Metoda Momentów

Metoda momentów polega na dopasowaniu teoretycznych momentów rozkładu do momentów empirycznych (obliczonych na podstawie danych), poprzez odpowiednie wyznaczenie parametrów rozkładu.

Aby użyć ją w R wystarczy użyć tej funkcji fitdist() z argumentem method="mme".

wek=rpois(100, 5)
fitdist(wek, "pois", method="mme")
## Fitting of the distribution ' pois ' by matching moments 
## Parameters:
##        estimate Std. Error
## lambda     4.93    2.22036

3.3 Ocena dopasowania rozkładu

Ocenę dopasowania można dokonać za pomocą:

  • testów dopasowania (o tym w następnym laboratorium),

  • kryteriów informacyjnych i statystyk dopasowania,

  • wykresów QQ-plot.

3.3.1 Kryteria informacyjne i statystyki dopasowania

Kryteria informacyjne, takie jak AIC (Akaike Information Criterion) i BIC (Bayesian Information Criterion), są używane do oceny jakości dopasowania modeli statystycznych, uwzględniając złożoność modelu. Pomagają one w wyborze najlepszego modelu spośród konkurencyjnych modeli, karząc za nadmierne dopasowanie (overfitting).

Niższe wartości AIC i BIC wskazują na lepsze dopasowanie modelu przy uwzględnieniu złożoności modelu.

Statystyki dopasowania są miarami używanymi do oceny, jak dobrze teoretyczny rozkład pasuje do danych empirycznych. Wyróżniamy między innymi:

  • Statystyka Kołmogorowa-Smirnowa mierzy maksymalną różnicę między dystrybuantą teoretyczną a empiryczną. Niższa wartość oznacza lepsze dopasowanie.

  • Statystyka Craméra-von Misesa mierzy sumę odległości między dystrybuantą teoretyczną a empiryczną. Niższa wartość oznacza lepsze dopasowanie.

  • Statystyka Andersona-Darlinga mierzy dopasowanie z większym naciskiem na krańce rozkładu. Niższe wartości wskazują na lepsze dopasowanie.

Funkcja gofstat() z biblioteki fitdistrplus dostarcza zestawienie różnych statystyk dopasowania dla każdego z dopasowanych rozkładów. Wystarczy, że jako argument podamy obiekt klasy fitdist lub listę obiektów fitdist.

Rozważmy poniższy przykład:

library(fitdistrplus)
set.seed(324)

Y_dane <- rgamma(100, 1,2)

fit1 <- fitdist(Y_dane, "norm")  # Dopasowanie rozkładu normalnego
fit2 <- fitdist(Y_dane, "exp")   # Dopasowanie rozkładu wykładniczego

# Porównanie jakości dopasowania obu rozkładów
gofstat(list(fit1, fit2),  fitnames = c("Norm", "Exp"))
## Goodness-of-fit statistics
##                                   Norm        Exp
## Kolmogorov-Smirnov statistic 0.1410153 0.07943591
## Cramer-von Mises statistic   0.5928108 0.04883928
## Anderson-Darling statistic   3.6702547 0.28922675
## 
## Goodness-of-fit criteria
##                                    Norm      Exp
## Akaike's Information Criterion 124.5233 56.14095
## Bayesian Information Criterion 129.7337 58.74612

W rezultacie otrzymujemy, że rozkład wykładniczy znacznie lepiej opisuje dane niż rozkład normalny, zarówno pod względem statystyk dopasowania, jak i kryteriów informacyjnych, ponieważ dla wszystkich przyjmował mniejsze wartości.

3.3.2 Wykresy QQ plot

QQ-plot (Quantile-Quantile plot) jest narzędziem do wizualnej oceny, jak dobrze dane empiryczne dopasowują się do wybranego teoretycznego rozkładu. Wykres ten porównuje kwantyle empiryczne danych z kwantylami teoretycznymi.

Możemy to zrobić ręcznie postępując według poniższych kroków:

  1. Posortuj zbiór danych w kolejności rosnącej.

  2. Wybierz rozkład wraz z parametrami, który chcesz sprawdzić, czy pasuje do posiadanego zbioru danych. Można użyć na przykład funkcji fitdist().

  3. Oblicz kwantyle teoretyczne dla wybranego rozkładu.

  4. Następnie narysuj wykres, dla którego na osi poziomej będą teoretyczne wartości kwantyli, a na osi pionowej posortowane obserwacje z kroku 1.

Jeśli punkty na wykresie znajdują się w przybliżeniu wzdłuż linii prostej, sugeruje to, że zestaw danych jest zgodny z założonym rozkładem. Odchylenia od linii prostej wskazują na odstępstwa od założonego rozkładu i wymagają dalszych badań.

y <- rgamma(200, shape = 2, rate = 0.5)
y_sorted <- sort(y)

# parametry dla rozkładu log-normal
y_fit=fitdist(y, distr = "lnorm")
# Obliczenie teoretycznych kwantyli dla rozkładu log normal
p <- ppoints(length(y)) # Generuje równomiernie rozłożone wartości kwantyli
theoretical_quantiles <- qlnorm(p, meanlog = y_fit$estimate[1], sdlog = y_fit$estimate[2])

plot( theoretical_quantiles, y_sorted)
abline(a=0,b=1,col="red")

Możemy również użyć funkcji qqcomp() z biblioteki fitdistrplus, która działa na obiektach klasy fitdist. Jako argument przekazujemy obiekt klasy fitdist lub listę obiektów klasy fitdist. Jeśli punkty na QQ plot układają się wzdłuż linii prostej \(y=x\), możemy uznać, że dopasowanie jest dobre. Jeśli nie, może to sugerować, że wybrany rozkład nie jest odpowiedni dla naszych danych.

y <- rgamma(200, shape = 2, rate = 0.5)
y_fit=fitdist(y, distr = "lnorm")
qqcomp(y_fit) # widzimy, że zdecydowanie pod koniec obserwacje nie leżą na prostej, zatem nie jest to prawdopodobnie dobry rozkład.

4 Przykład - Analiza danych

Ładujemy zbiór danych faithful z biblioteki fitdistrplus. Zbiór danych zapiera informacje o czasie trwania erupcji (w minutach) dla gejzeru Old Faithful w Parku Narodowym Yellowstone, Wyoming, USA oraz czasami oczekiwań (w minutach) między erupcjami. Wykonaj poniższy kod.

library(fitdistrplus)
attach(faithful)
summary(faithful)
value <- eruptions
hist(value, breaks = 30)
boxplot(value) # Zauważmy, że dużo wartości odstających. Zatem, prawdopodobnie będzie to musiał być rozkład który jest zdefiniowany tylko dla wartości dodatnich i ma `ciężki` ogon.

descdist(value, discrete = FALSE) # Funkcja pomagająca w dopasowaniu rozkładu do danych

normalize <- function(x) {
  return((x - min(x)) / (max(x) - min(x)))
}

# Rozkład beta jest określony na przedziale otwartym (0,1) więc musimy wykonać normalizacje danych i przesunąć delikatnie skrajne wartości tak aby był to przedział otwarty.
normalized_value <- normalize(value)
normalized_value[normalized_value==0]=0.00001
normalized_value[normalized_value==1]=0.99999

fit_beta <- fitdist(normalized_value, "beta")
plot(fit_beta) # Wyświetla wykres dopasowania 

fit_unif <- fitdist(normalized_value, "unif")
plot(fit_unif) # Wyświetla wykres dopasowania 


#Ocena dopasowania
gof_results <- gofstat(list(fit_beta, fit_unif), fitnames = c("Beta", "Jedostajny"))
print(gof_results)

Wnioski: Rozkład beta jest bardziej odpowiedni dla danych na podstawie kryteriów informacyjnych, ponieważ AIC, jak i BIC są niższe dla rozkładu beta, niż dla rozkładu jednostajnego.

Dopasowanie rozkładu do prawdziwych danych jest często wyzwaniem z wielu powodów. Dane empiryczne mogą mieć złożone struktury, takie jak wielomodalność (może oznaczać, że dane mogą pochodzić z mieszaniny kilku różnych rozkładów), asymetrię (która często utrudnia dopasowanie rozkładów symetrycznych, takich jak normalny), czy obecność wartości odstających (mogące znacząco wpływać na dopasowanie modelu, prowadząc do nieadekwatnych parametrów), które utrudniają dopasowanie prostych modeli teoretycznych.

Zatem można rozpatrywać:

4.1 Zadanie - analiza danych

Wykonaj podobną analizę danych jak w powyższym przykładzie dla wybranych danych ilościowych ze zbioru danych daneSoc z biblioteki Przewodnik. Napisz na koniec wniosek z wykonanej analizy.

Podpowiedź: zaczniemy od danych, dla każdej zmiennej ds, dr, dw przeprowadzić osobną analizę:

library("Przewodnik")
dane <- daneSoc
summary(dane)
ds <- dane$cisnienie.skurczowe #wypróbować rozkład normalny `norm`
dr <- dane$cisnienie.rozkurczowe #wypróbować rozkład logistyczni `logis` oraz normalny `norm`
dw <- dane$wiek #wypróbować rozkład beta `beta` oraz jednostajny `unif`