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
symulowanie losowości, 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 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 łatwa implementacja 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.
Algorytm generuje zatem liczby z przedziału \([0,m-1]\).
1.2 Zadanie 1 - generator
liczb losowych
Zaimplementuj funkcję lcm(N, 0, 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?
1.3 Zadanie 2 - generator
liczb losowych
Wywołaj wcześniej napisana funkcję dla poniższych argumentów. Czy w
tym przypadku otrzymane wartości są bardziej losowe?
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:
Możliwe wartości dla prefixów to:
r (random) - generuje próbę prostą z danego rozkładu
o rozmiarze n (pierwszy argument funkcji),
p (probability) - wyznacza wartość dystrybuanty w
punktach określonych przez wektor x (pierwszy argument
funkcji),
d (density) - wyznacza wartość gęstości (dla
rozkładów ciągłych) lub prawdopodobieństwa (dla rozkładów dyskretnych)
danego rozkładu w punktach określonych przez wektor x
(pierwszy argument funkcji),
q (quantile) - wyznacza kwantyl danego rozkładu w
punktach x (pierwszy argument funkcji).
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:
## [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:
2.1.1 Zadania - rozkład
jednostajny
2.1.1.1 Zadanie 1
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.
Przykładowy wykres
2.1.1.2 Zadanie 2
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\).
Przykładowy wykres
2.1.1.3 Zadanie 3
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]\).
Przykładowy wykres
2.1.2 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:
znana,
ciągła,
ś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]\)).
Dowód:
Niech \(F_X\) oznacza dystrybuantę
zmiennej losowej \(X=F^{-1}(U)\).
Wtedy, \[F_X(x)=P(X\leq
x)=P\left(F^{-1}(U)\leq x\right)=P\left(F\left(F^{-1}(U)\right)\leq
F(x)\right)=P\left(U\leq F(x)\right)=F(x).\]
W rezultacie mamy następujący algorytm na generowanie zmiennych
losowych \(X\) o rozkładzie \(F_X\):
wygeneruj niezależny zmienne losowe \(u_1, \dots, u_n\) z rozkładu jednostajnego
\(U[0,1]\),
oblicz \(x_1=F_X^{-1}(u_1), \dots,
x_n=F_X^{-1}(u_n)\).
Wtedy zmienne losowe \(x_1, \dots,
x_n\) są niezależne o rozkładzie \(F_X\).
2.1.3 Zadanie 1 - metoda
dystrybuanty odwrotnej
Wygeneruj 5 obserwacji z rozkłady 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.\]
2.1.4 Zadanie 2 - metoda
dystrybuanty odwrotnej
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.
2.1.5 Uwagi o metodzie
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.2 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")
2.2.1 Zadania - rozkład
normaly
2.2.1.1 Zadanie 1
Przetestuj powyższy kod dla innych parametrów rozkładu normalnego.
Jak zmiana parametrów wpływa na kształt wykresów?
2.2.1.2 Zadanie 2
Załóżmy, że wyniki egzaminu przeprowadzonego dla \(200\) studentów, mają rozkład normalny o
średniej \(72\) i odchyleniu
standardowym \(15.2\). Wygeneruj wyniki
z egzaminu i wyznacz jakie jest prawdopodobieństwo uzyskania \(84\) punkty lub więcej na egzaminie.
2.2.1.3 Zadanie 3
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)\)).
2.3 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.3.1 Zadania - 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.
2.4 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.4.1 Zadania - rozkład
t-Studenta
2.4.1.1 Zadanie 1
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.4.1.2 Zadanie 2
Wygeneruj \(10\) obserwacji z
rozkładu normalnego o średniej \(0.5\)
i wariancji \(2\), a następnie oblicz
statystykę \(t\). Narysuj rozkład
t-Studenta o \(9\) stopniach swobody i
zaznacz obliczoną wartość statystyki \(t\) na wykresie.
2.5 Rozkład F
Rozkład F nazywany jest również rozkładem F - Snedecora lub rozkładem
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 wartości 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")
2.6 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 R nazywa się
binom.
Wygenerowanie zmiennych z rozkładu Bernoulliego przy pomocy rozkładu
jednostajnego \(U[0,1]\) jest również
możliwe.
Algorytm:
Generujemy liczby \(u_1,\dots,
u_n\) z rozkładu \(U[0,1]\).
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.6.1 Zadania - rozkład
Bernoullego
2.6.1.1 Zadanie 1
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.
2.6.1.2 Zadanie 2
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 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 wektora \(x\) do 500 i
powtórz symulację.
2.7 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 zdarzeń zaobserwowanych przez detektor w ciągu
minuty,
liczba włączeń klimatyzacji w ciągu godziny,
liczba budynków na kilometr kwadratowy,
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.
2.7.1 Zadanie - rozkład
Poissona
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.
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:
opisaniu danych - poprawne zrozumienie i opisanie właściwości
danych,
przewidywania,
testowania hipotez,
symulacji - tworzenie realistycznych symulacji na podstawie
modelu.
Dla danego zestawu danych powinno się wykonać następujące kroki:
Znalezienie kandydata na rozkład.
Estymacja parametrów rozkładu.
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 z 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 (mean) i błąd standardowy dla tych parametrów
(sd).
## Fitting of the distribution ' norm ' by maximum likelihood
## Parameters:
## estimate Std. Error
## mean 0.2353006 0.2977325
## sd 2.9773248 0.2105285
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".
## Fitting of the distribution ' pois ' by matching moments
## Parameters:
## estimate Std. Error
## lambda 5.18 0.2275961
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:
Posortuj zbiór danych w kolejności rosnącej.
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().
Oblicz kwantyle teoretyczne dla wybranego rozkładu.
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ć:
Mieszanki rozkładów (mixture models), czyli kombinacje kilku
rozkładów, które razem modelują dane.
W niektórych przypadkach, transformacja danych może uprościć
dopasowanie. Przykłady to transformacja Box-Cox lub
log-transformacja.
Po dopasowaniu modelu, analiza reszt może pomóc zidentyfikować
problemy z dopasowaniem i wskazać kierunki do dalszej poprawy
modelu.
W niektórych przypadkach przydatne może się okazać zastosowanie
metod nieliniowych, takich jak krzywe nieliniowe lub regresja
nieliniowa.
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.
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 symulowanie losowości, 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 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 łatwa implementacja 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.
Algorytm generuje zatem liczby z przedziału \([0,m-1]\).
1.2 Zadanie 1 - generator liczb losowych
Zaimplementuj funkcję lcm(N, 0, 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?
1.3 Zadanie 2 - generator liczb losowych
Wywołaj wcześniej napisana funkcję dla poniższych argumentów. Czy w tym przypadku otrzymane wartości są bardziej losowe?
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:
Możliwe wartości dla prefixów to:
r(random) - generuje próbę prostą z danego rozkładu o rozmiarzen(pierwszy argument funkcji),p(probability) - wyznacza wartość dystrybuanty w punktach określonych przez wektorx(pierwszy argument funkcji),d(density) - wyznacza wartość gęstości (dla rozkładów ciągłych) lub prawdopodobieństwa (dla rozkładów dyskretnych) danego rozkładu w punktach określonych przez wektorx(pierwszy argument funkcji),q(quantile) - wyznacza kwantyl danego rozkładu w punktachx(pierwszy argument funkcji).
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:
## [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:
2.1.1 Zadania - rozkład
jednostajny
2.1.1.1 Zadanie 1
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.
Przykładowy wykres
2.1.1.2 Zadanie 2
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\).
Przykładowy wykres
2.1.1.3 Zadanie 3
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]\).
Przykładowy wykres
2.1.2 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:
znana,
ciągła,
ś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]\)).
Dowód:
Niech \(F_X\) oznacza dystrybuantę
zmiennej losowej \(X=F^{-1}(U)\).
Wtedy, \[F_X(x)=P(X\leq
x)=P\left(F^{-1}(U)\leq x\right)=P\left(F\left(F^{-1}(U)\right)\leq
F(x)\right)=P\left(U\leq F(x)\right)=F(x).\]
W rezultacie mamy następujący algorytm na generowanie zmiennych
losowych \(X\) o rozkładzie \(F_X\):
wygeneruj niezależny zmienne losowe \(u_1, \dots, u_n\) z rozkładu jednostajnego
\(U[0,1]\),
oblicz \(x_1=F_X^{-1}(u_1), \dots,
x_n=F_X^{-1}(u_n)\).
Wtedy zmienne losowe \(x_1, \dots,
x_n\) są niezależne o rozkładzie \(F_X\).
2.1.3 Zadanie 1 - metoda
dystrybuanty odwrotnej
Wygeneruj 5 obserwacji z rozkłady 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.\]
2.1.4 Zadanie 2 - metoda
dystrybuanty odwrotnej
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.
2.1.5 Uwagi o metodzie
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.2 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")
2.2.1 Zadania - rozkład
normaly
2.2.1.1 Zadanie 1
Przetestuj powyższy kod dla innych parametrów rozkładu normalnego.
Jak zmiana parametrów wpływa na kształt wykresów?
2.2.1.2 Zadanie 2
Załóżmy, że wyniki egzaminu przeprowadzonego dla \(200\) studentów, mają rozkład normalny o
średniej \(72\) i odchyleniu
standardowym \(15.2\). Wygeneruj wyniki
z egzaminu i wyznacz jakie jest prawdopodobieństwo uzyskania \(84\) punkty lub więcej na egzaminie.
2.2.1.3 Zadanie 3
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)\)).
2.3 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.3.1 Zadania - 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.
2.4 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.4.1 Zadania - rozkład
t-Studenta
2.4.1.1 Zadanie 1
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.4.1.2 Zadanie 2
Wygeneruj \(10\) obserwacji z
rozkładu normalnego o średniej \(0.5\)
i wariancji \(2\), a następnie oblicz
statystykę \(t\). Narysuj rozkład
t-Studenta o \(9\) stopniach swobody i
zaznacz obliczoną wartość statystyki \(t\) na wykresie.
2.5 Rozkład F
Rozkład F nazywany jest również rozkładem F - Snedecora lub rozkładem
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 wartości 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")
2.6 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 R nazywa się
binom.
Wygenerowanie zmiennych z rozkładu Bernoulliego przy pomocy rozkładu
jednostajnego \(U[0,1]\) jest również
możliwe.
Algorytm:
Generujemy liczby \(u_1,\dots,
u_n\) z rozkładu \(U[0,1]\).
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.6.1 Zadania - rozkład
Bernoullego
2.6.1.1 Zadanie 1
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.
2.6.1.2 Zadanie 2
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 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 wektora \(x\) do 500 i
powtórz symulację.
2.7 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 zdarzeń zaobserwowanych przez detektor w ciągu
minuty,
liczba włączeń klimatyzacji w ciągu godziny,
liczba budynków na kilometr kwadratowy,
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.
2.7.1 Zadanie - rozkład
Poissona
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.
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:
## [1] 0.2002145 0.6852186 0.9168758 0.2843995 0.1046501 0.7010575 0.5279600
## [8] 0.8079352 0.9565001 0.1104530Funkcja 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:
2.1.1 Zadania - rozkład jednostajny
2.1.1.1 Zadanie 1
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.
Przykładowy wykres
2.1.1.2 Zadanie 2
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\).
Przykładowy wykres
2.1.1.3 Zadanie 3
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]\).
Przykładowy wykres
2.1.2 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:
znana,
ciągła,
ś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]\)).
Dowód:
Niech \(F_X\) oznacza dystrybuantę zmiennej losowej \(X=F^{-1}(U)\). Wtedy, \[F_X(x)=P(X\leq x)=P\left(F^{-1}(U)\leq x\right)=P\left(F\left(F^{-1}(U)\right)\leq F(x)\right)=P\left(U\leq F(x)\right)=F(x).\]
W rezultacie mamy następujący algorytm na generowanie zmiennych losowych \(X\) o rozkładzie \(F_X\):
wygeneruj niezależny zmienne losowe \(u_1, \dots, u_n\) z rozkładu jednostajnego \(U[0,1]\),
oblicz \(x_1=F_X^{-1}(u_1), \dots, x_n=F_X^{-1}(u_n)\).
Wtedy zmienne losowe \(x_1, \dots, x_n\) są niezależne o rozkładzie \(F_X\).
2.1.3 Zadanie 1 - metoda dystrybuanty odwrotnej
Wygeneruj 5 obserwacji z rozkłady 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.\]
2.1.4 Zadanie 2 - metoda dystrybuanty odwrotnej
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.
2.1.5 Uwagi o metodzie 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.2 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")
2.2.1 Zadania - rozkład normaly
2.2.1.1 Zadanie 1
Przetestuj powyższy kod dla innych parametrów rozkładu normalnego. Jak zmiana parametrów wpływa na kształt wykresów?
2.2.1.2 Zadanie 2
Załóżmy, że wyniki egzaminu przeprowadzonego dla \(200\) studentów, mają rozkład normalny o średniej \(72\) i odchyleniu standardowym \(15.2\). Wygeneruj wyniki z egzaminu i wyznacz jakie jest prawdopodobieństwo uzyskania \(84\) punkty lub więcej na egzaminie.
2.2.1.3 Zadanie 3
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)\)).
2.3 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.3.1 Zadania - 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.
2.4 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.4.1 Zadania - rozkład t-Studenta
2.4.1.1 Zadanie 1
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.4.1.2 Zadanie 2
Wygeneruj \(10\) obserwacji z rozkładu normalnego o średniej \(0.5\) i wariancji \(2\), a następnie oblicz statystykę \(t\). Narysuj rozkład t-Studenta o \(9\) stopniach swobody i zaznacz obliczoną wartość statystyki \(t\) na wykresie.
2.5 Rozkład F
Rozkład F nazywany jest również rozkładem F - Snedecora lub rozkładem 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 wartości 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")
2.6 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 R nazywa się
binom.
Wygenerowanie zmiennych z rozkładu Bernoulliego przy pomocy rozkładu jednostajnego \(U[0,1]\) jest również możliwe.
Algorytm:
Generujemy liczby \(u_1,\dots, u_n\) z rozkładu \(U[0,1]\).
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.6.1 Zadania - rozkład Bernoullego
2.6.1.1 Zadanie 1
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.
2.6.1.2 Zadanie 2
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 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 wektora \(x\) do 500 i
powtórz symulację.
2.7 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 zdarzeń zaobserwowanych przez detektor w ciągu minuty,
liczba włączeń klimatyzacji w ciągu godziny,
liczba budynków na kilometr kwadratowy,
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.
2.7.1 Zadanie - rozkład Poissona
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.
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:
opisaniu danych - poprawne zrozumienie i opisanie właściwości danych,
przewidywania,
testowania hipotez,
symulacji - tworzenie realistycznych symulacji na podstawie modelu.
Dla danego zestawu danych powinno się wykonać następujące kroki:
Znalezienie kandydata na rozkład.
Estymacja parametrów rozkładu.
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 z 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 (mean) i błąd standardowy dla tych parametrów
(sd).
## Fitting of the distribution ' norm ' by maximum likelihood
## Parameters:
## estimate Std. Error
## mean 0.2353006 0.2977325
## sd 2.9773248 0.2105285
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".
## Fitting of the distribution ' pois ' by matching moments
## Parameters:
## estimate Std. Error
## lambda 5.18 0.2275961
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:
Posortuj zbiór danych w kolejności rosnącej.
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().
Oblicz kwantyle teoretyczne dla wybranego rozkładu.
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.
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 z 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 (mean) i błąd standardowy dla tych parametrów
(sd).
## Fitting of the distribution ' norm ' by maximum likelihood
## Parameters:
## estimate Std. Error
## mean 0.2353006 0.2977325
## sd 2.9773248 0.21052853.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".
## Fitting of the distribution ' pois ' by matching moments
## Parameters:
## estimate Std. Error
## lambda 5.18 0.22759613.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.74612W 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:
Posortuj zbiór danych w kolejności rosnącej.
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().Oblicz kwantyle teoretyczne dla wybranego rozkładu.
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ć:
Mieszanki rozkładów (mixture models), czyli kombinacje kilku rozkładów, które razem modelują dane.
W niektórych przypadkach, transformacja danych może uprościć dopasowanie. Przykłady to transformacja Box-Cox lub log-transformacja.
Po dopasowaniu modelu, analiza reszt może pomóc zidentyfikować problemy z dopasowaniem i wskazać kierunki do dalszej poprawy modelu.
W niektórych przypadkach przydatne może się okazać zastosowanie metod nieliniowych, takich jak krzywe nieliniowe lub regresja nieliniowa.
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.