16 KiB
Modelowanie przestrzeni probabilistycznych w języku R cz. 2
Kontynuujemy temat generowania przestrzeni probabilistycznych w języku R. Przypominamy, że większość interesujących nas funkcji znajduje się w pakiecie probs.
Prawdopodobieństwo warunkowe
Przypomnijmy, że jeśli dane są dwa zdarzenia $A$ i $B$ określone w tej samej przestrzeni probabilistycznej oraz $\mathbb{P}(B)>0$, to możemy zdefiniować prawdopodobieństwo warunkowe zdarzenia $A$ pod warunkiem zajścia zdarzenia $B$ w następujący sposób:
$$\mathbb{P}(A|B)=\frac{\mathbb{P}(A\cap B)}{\mathbb{P}(B)}.$$
Jeśli chcemy obliczyć takie prawdopodobieństwo w języku R, możemy użyć polecenia
Prob(x, event=A, given=B),
gdzie w miejscu $A$ i $B$ należy podać odpowiednie podzbiory w formie ramek danych lub formuły logicznej opisującej to zdarzenie. Jeśli $A$ zostanie wcześniej zdefiniowane jako podzbiór odpowiedniej przestrzeni probabilistycznej, to można też użyć polecenia Prob(A, given=B).
Przykład 1
Rzucamy dwukrotnie sześcienną kostką do gry. Ile wynosi prawdopodobieństwo, że suma wyrzuconych oczek jest większa od $7$, jeśli wiemy, że w pierwszym rzucie wypadła parzysta liczba oczek?
# Definiujemy przestrzeń probabilistyczną odpowiadającą temu doświadczeniu losowemu
kostki2 = rolldie(2, makespace=TRUE)
# Definiujemy zdarzenia A i B. Skorzystamy z operatora %% (x %% y zwraca resztę z dzielenia x przez y)
A = subset(kostki2, X1 %% 2 == 0)
B = subset(kostki2, X1 + X2 > 7)
# Obliczamy prawdopodobieństwo warunkowe P(A|B)
p = Prob(A, given=B)
print(p)
# Możemy porównać to prawdopodobieństwo do prawdopodobieństwa zdarzenia A
p2 = Prob(A)
print(p2)Przestrzenie produktowe
Przypomnijmy, że jeśli dane są dwie przestrzenie probabilistyczne $(\Omega_1, \mathcal{F}_1, \mathbb{P}_1)$ oraz $(\Omega_2, \mathcal{F}_2, \mathbb{P}_2)$, to możemy zdefiniować ich produkt jako przestrzeń probabilistyczną $(\Omega,\mathcal{F},\mathcal{P})$, gdzie $\Omega = \Omega_1\times\Omega_2$, $\mathcal{F}$ zawiera wszystkie zbiory postaci $A_1\times A_2\in \mathcal{F}_1\times\mathcal{F}_2$, oraz funkcja prawdopodobieństwa $\mathbb{P}$ spełnia warunek
$$ \mathbb{P}(A_1\times A_2) = \mathbb{P}_1(A_1)\cdot\mathbb{P}_2(A_2).$$
W analogiczny sposób możemy definiować produkt dowolnej, skończonej liczby przestrzeni probabilistycznych. W języku R istnieje kilka sposobów modelowania tego typu przestrzeni.
Niezależne powtórzenia tego samego eksperymentu losowego
Załóżmy najpierw, że wielokrotnie w sposób niezależny powtarzamy ten sam eksperyment losowy, np. rzucamy $n$-krotnie standardową kostką sześcienną. Wówczas możemy użyć polecenia
iidspace(x, ntrials, probs=NULL),
gdzie x oznacza wektor możliwych wyników w pojedynczym powtórzeniu eksperymentu losowego, ntrials podaje liczbę powtórzeń, a probs to wektor prawopodobieństw odpowiadający wynikom x. Jeśli wszystkie wyniki pojedynczego powtórzenia eksperymentu losowego są równo prawdopodobne, to możemy nie podawać wektora probs.
Przykład 2
Wygeneruj przestrzeń probabilistyczną odpowiadającą trzykrotnemu rzutowi czworościenną kostką przy założeniu, że:
- kostka jest symetryczna (tzn. każdy wynik jest równo prawdopodobny),
- prawdopodobieństwo wypadnięcia $1$ to $1/2$, a prawdopodobieństwo wypadnięcia dowolnej innej liczby oczek to $1/6$.
# Generujemy przestrzeń probabilistyczną odpowiadającą 3-krotnemu rzutowi symetryczną kostką
x = 1:4
kostka_sym = iidspace(x, 3)
print(kostka_sym)
# Zwróćmy uwagę, że równoważnie można było w tym przypadku użyć polecenia rolldie(3, nsides=4, makespace=TRUE)
# Generujemy przestrzeń probabilistyczną odpowiadającą 3-krotnemu rzutowi ,,oszukaną'' kostką
kostka_osz = iidspace(x, 3, probs = c(1/2, 1/6, 1/6, 1/6))
print(kostka_osz)Szczególnym przypadkiem niezależnych powtórzeń tego samego eksperymentu losowego jest sytuacja, gdy ten eksperyment losowy ma tylko dwa możliwe wyniki umownie nazywane ,,sukcesem'' i ,,porażką''. Jeśli liczba powtórzeń eksperymentu losowego jest ustalona i skończona, to mamy wówczas do czynienia z tzw. schematem Bernoulliego. Prawdopodobieństwo uzyskania dokładnie $k$ sukcesów w tego typu doświadczeniu losowym jest wtedy równe:
$$\tau_k = {n\choose k}p^k(1-p)^{n-k} \ \text{ dla } \ k=0,1,\ldots,n,$$
gdzie $n$ to liczba powtórzeń, a $p$ to prawdopodobieństwo sukcesu w pojedynczej próbie. W R prawdopodobieństwa tego typu można liczyć przy pomocy następującego polecenia z pakietu stats:
dbinom(x, size, prob),
gdzie x oznacza wektor zawierający interesujące nas liczby sukcesów, size oznacza liczbę powtórzeń, a prob prawdopodobieństwo sukcesu w pojedynczej próbie.
Przykład 3
Losujemy $10$ razy, ze zwracaniem jedną kartę z talii $52$ kart. Ile wynosi prawdopodobieństwo, że wylosujemy:
- dokładnie $6$ kierów,
- przynajmniej $6$ kierów.
Uwaga: Poniższy przykład wymaga zainstalowania i załadowania pakietu stats.
# Zwróćmy uwagę, że ponieważ losujemy karty ze zwracaniem, to mamy do czynienia z niezależnymi powtórzeniami tego samego eksperymentu losowego
# Interesuje nas liczba kierów, więc jako sukces możemy przyjąć wylosowanie kiera
# Liczba powtórzeń wynosi n=10, a prawdopodobieństwo sukcesu w pojedynczej próbie to p=1/4 (bo w talii mamy 13 kierów)
p1 = dbinom(6, 10, 1/4)
print(paste('Prawdopodobieństwo wylosowania dokładnie 6 kierów: ', p1))
# W przypadku zdarzenia ,,wylosujemy przynajmniej 6 kierów'' jako x podamy wektor zawierający wszystkie interesujące nas liczby sukcesów
y = dbinom(6:10, 10, 1/4)
print(y)
# Zwróćmy uwagę, że zostanie zwrócony wektor prawdopodobieństw odpowiadający wszystkim możliwym liczbom sukcesów. Aby uzyskać odpowiedź musimy je zsumować
p2 = sum(y)
print(paste('Prawdopodobieństwo wylosowania przynajmniej 6 kierów: ', p2))[1] "Prawdopodobieństwo wylosowania dokładnie 6 kierów: 0.0162220001220703" [1] 1.622200e-02 3.089905e-03 3.862381e-04 2.861023e-05 9.536743e-07 [1] "Prawdopodobieństwo wylosowania przynajmniej 6 kierów: 0.0197277069091797"
Mogą nas też interesować doświadczenia losowe polegające na powtarzaniu eksperymentu losowego o dwóch wynikach (ponownie umownie nazywanymi "sukcesem" i "porażką") do momentu osiągnięcia pierwszego sukcesu. Mamy wówczas do czynienia z tzw. rozkładem geometrycznym i możemy obliczyć prawdopodobieństwo, że będziemy potrzebowali dokładnie $k$ prób na osiągnięcie pierwszego sukcesu zgodnie ze wzorem:
$$\sigma_k=p(1-p)^{k-1} \ \text{ dla } \ k=1,2,\ldots,$$
gdzie $p$ oznacza prawdopodobieństwo sukcesu w pojedynczej próbie. W R prawdopodobieństwa tego typu można liczyć przy pomocy następującego polecenia z pakietu stats:
dgeom(x, prob),
gdzie x oznacza wektor zawierający interesujące nas liczby porażek przed osiągnięciem pierwszego sukcesu, a prob prawdopodobieństwo sukcesu w pojedynczej próbie.
Przykład 4
Losujemy wielokrotnie ze zwracaniem jedną kartę z talii $52$ kart do momentu wylosowania pierwszego kiera. Ile wynosi prawdopodobieństwo, że wykonamy:
- dokładnie $6$ losowań,
- co najwyżej $6$ losowań,
- co najmniej $7$ losowań.
Uwaga: Poniższy przykład wymaga zainstalowania i załadowania pakietu stats.
# Podobnie, jak w poprzednim przykładzie mamy do czynienia z niezależnymi powtórzeniami tego samego eksperymentu losowego, a sukcesem jest wylosowanie kiera
# Zwróćmy uwagę, że w poleceniu dgeom podajemy liczbę porażek przed pierwszym sukcesem
p1 = dgeom(5, 1/4)
print(paste('Prawdopodobieństwo dokładnie 6 losowań: ', p1))
# W przypadku zdarzenia ,,wykonamy co najwyżej 6 losowań'' jako x podamy wektor zawierający wszystkie interesujące nas liczby porażek
y = dgeom(0:5, 1/4)
print(y)
# Zwróćmy uwagę, że zostanie zwrócony wektor prawdopodobieństw odpowiadający wszystkim możliwym liczbom wykonanych losowań. Aby uzyskać odpowiedź musimy je zsumować
p2 = sum(y)
print(paste('Prawdopodobieństwo co najwyżej 6 losowań: ', p2))
# W przypadku zdarzenia ,,wykonamy co najmniej 7 losowań'' nie będziemy tego rozbijać na wszystkie możliwe przypadki (bo byłoby ich nieskończenie wiele). Zamiast tego skorzystamy ze zdarzenia przeciwnego
p3 = 1-p2
print(paste('Prawdopodobieństwo co najmniej 7 losowań: ', p3))[1] "Prawdopodobieństwo dokładnie 6 losowań: 0.059326171875" [1] 0.25000000 0.18750000 0.14062500 0.10546875 0.07910156 0.05932617 [1] "Prawdopodobieństwo co najwyżej 6 losowań: 0.822021484375" [1] "Prawdopodobieństwo co najmniej 7 losowań: 0.177978515625"
Ogólne przestrzenie produktowe
Na koniec pokażemy, jak możemy modelować przestrzenie produktowe w przypadku, gdy łączymy wyniki kilku różnych eksperymentów losowych. Wiemy już, że w takim przypadku zbiór zdarzeń elementarnych $\Omega$ jest iloczynem kartezjańskim zbiorów zdarzeń elementarnych odpowiadających poszczególnym eksperymentom losowym. W języku R do generowania iloczynu kartezjańskiego wektorów możemy użyć polecenia expand.grid. Jako argumenty należy podać wektory, których iloczyn kartezjański chcemy wyznaczyć. Jeśli podamy je w formie tag=v, to kolumna odpowiadająca wektorowi v będzie miała nazwę tag.
Przykład 5
Rzucamy jednokrotnie sześcienną kostką i symetryczną monetą. Wygeneruj przestrzeń probabilistyczną odpowiadającą temu doświadczeniu losowemu.
# Najpierw generujemy zbiór zdarzeń elementarnych
Omega = expand.grid(k = 1:6, m = c('O', 'R'))
# Następnie na tej podstawie generujemy przestrzeń probabilistyczną. W tym przypadku możemy przyjąć, że mamy do czynienia z modelem klasycznym
X = probspace(Omega)Error in probspace(Omega): could not find function "probspace" Traceback:
Należy zwrócić uwagę, że jeśli definiujemy w R produktową przestrzeń probabilistyczną opierając się na przestrzeniach probabilistycznych, w których prawdopodobieństwo nie jest definiowane w sposób klasyczny, to musimy dodatkowo zdefiniować w odpowiedni sposób wektor prawdopodobieństw.
Przykład 6
Rzucamy jednokrotnie kostką i monetą. Zakładamy, że na kostce jedno oczko wypada z prawdopodobieństwem $1/2$, a pozostałe liczby oczek wypadają z prawdopodobieństwem $1/10$. Wygeneruj przestrzeń probabilistyczną odpowiadającą temu doświadczeniu losowemu. Następnie oblicz prawdopodobieństwo zdarzenia $A$ - wypadł orzeł i nieparzysta liczba oczek.
# Możemy wykorzystać zbiór zdarzeń elementarnych Omega wygenerowany w poprzednim przykładzie
# Zgodnie z definicją przestrzeni produktowej aby znaleźć prawdopodobieństwo zdarzenia elementarnego (x,y) należy pomnożyć prawdopodobieństwa tych zdarzeń w pierwotnych przestrzeniach probabilistycznych
p = rep(0, 12)
for (i in 1:12) {
if (Omega[i,1] == 1) {p[i] = 0.5*0.5} else {p[i] = 0.1*0.5}
}
# Generujemy przestrzeń probabilistyczną
X2 = probspace(Omega, p)
print(X2)
# Możemy teraz obliczyć prawdopodobieństwo zdarzenie A
A = subset(Omega, (k %in% seq(1, 6, by=2)) & (m == 'O'))
print(Prob(A))
# Alternatywnie, można by wyznaczyć prawdopodobieństwo wypadnięcia nieparzystej liczby oczek w przestrzeni rzutu kostką i prawdopodobieństwo wypadnięcia orła w przestrzeni rzutu monetą i pomnożyć przez siebie te dwa prawdopodobieństwa
kostka = probspace(1:6, c(0.5, rep(0.1,5))) # Definiujemy przestrzeń związaną z rzutem kostką
A1 = subset(kostka, (x %in% seq(1, 6, by=2)))
moneta = probspace(c('O', 'R')) # Definiujemy przestrzeń związaną z rzutem monetą
A2 = subset(moneta, x == 'O')
pA = Prob(A1) * Prob(A2)
print(pA)