sample(c("H", "T"), 10, replace = TRUE) [1] "H" "T" "T" "T" "H" "T" "T" "T" "T" "T"Scopul acestor note de laborator este de a introduce o serie de noțiuni și funcții utile în rezolvarea problemelor de probabilități prin intermediul programului R.
Modelarea problemelor de probabilități
În această secțiune vom prezenta o serie de exerciții și probleme de probabilități, de dificultăți diferite, pentru a ilustra cum pot fi ele rezolvate prin intermediul limbajului R.
Aruncatul cu banul (I)
În acest exemplu vrem să simulăm aruncarea unei monede (echilibrate) folosind funcția sample(). Această funcție permite extragerea, cu sau fără întoarcere (prin intermediul argumentului replace = TRUE sau replace = FALSE - aceasta este valoarea prestabilită), a unui eșantion de volum dat (size) dintr-o mulțime de elemente x.
Spre exemplu dacă vrem să simulăm \(10\) aruncări cu banul atunci apelăm:
Pentru a estima probabilitatea de apariției a stemei (H) repetăm aruncarea cu banul de \(10000\) de ori și calculăm raportul dintre numărul de apariții ale evenimentului \(A=\{H\}\) și numărul total de aruncări:
# atunci cand moneda este echilibrata
ban <- sample(c("H","T"), 10000, replace = TRUE)
table(ban)/length(ban)ban
H T
0.501 0.499 Pentru cazul în care moneda nu este echilibrată putem specifica prin argumentul prob =, probabilitatea de a obține cap respectiv coadă:
ban <- sample(c("H","T"), 10000, replace = TRUE, prob = c(0.2, 0.8))
table(ban)/length(ban)ban
H T
0.1949 0.8051 Putem vedea cum evoluează această frecvență relativă în funcție de numărul de aruncări
N <- 10000
ban <- sample(c("H","T"), N, replace = TRUE)
# sirul frecventelor relative
y <- cumsum(ban == "H")/(1:N)
plot(1:N, y, type = "o", col = myblue, bty = "n",
xlab ="Numar aruncari", ylab = "Frecventa relativa",
ylim = c(min(y), max(y)), pch = 16)
abline(h = 0.5, lty = 2, col = myred)
Cum puteți exprima în termeni de aruncat cu banul problema expusă la curs ?
Numărul de băieți dintr-o familie cu doi copii
În lucrarea Mathematical games din revista Scientific American numărul 200 (pag. 164-174), Martin Gardner a formulat următoarea problemă:
Exercițiul 1.1 O familie are doi copii. Care este probabilitatea ca ambii copii să fie băieți știind că cel puțin unul dintre copii este băiat ? Care este probabilitatea ca ambii copii să fie băieți știind că cel mai tânăr este băiat ?
Soluție. Pentru a răspunde la cele două întrebări să observăm că cei doi copii (copilul mai mare și cel mai mic) pot fi ambii de sex masculin sau feminin, prin urmare avem patru combinații de sexe, pe care le presupunem egal probabile. Putem reprezenta spațiul stărilor prin
\[ \Omega = \{BB,BF,FB,FF\} \]
unde \(\mathbb{P}(BB)=\mathbb{P}(BF)=\mathbb{P}(FB)=\mathbb{P}(FF)=\frac{1}{4}\).
Pentru a răspunde la prima întrebare avem:
\[ \begin{aligned} \mathbb{P}\left(BB\,|\,\text{cel puțin unul este băiat}\right) &= \mathbb{P}\left(BB\,|\,BF\cup FB\cup BB\right)\\ &= \frac{\mathbb{P}\left(BB\cap (BF\cup FB\cup BB)\right)}{\mathbb{P}\left(BF\cup FB\cup BB\right)} \\ &= \frac{\mathbb{P}\left(BB\right)}{\mathbb{P}\left(BF\cup FB\cup BB\right)} = \frac{1}{3}. \end{aligned} \]
Iar pentru cea de-a doua întrebare:
\[ \begin{aligned} \mathbb{P}\left(BB\,|\,\text{cel mai tânăr este băiat}\right) &= \mathbb{P}\left(BB\,|\,FB\cup BB\right)\\ &= \frac{\mathbb{P}\left(BB\cap (FB\cup BB)\right)}{\mathbb{P}\left(FB\cup BB\right)} \\ &= \frac{\mathbb{P}\left(BB\right)}{\mathbb{P}\left(FB\cup BB\right)} = \frac{1}{2}. \end{aligned} \]
Vom încerca să răspundem la aceste întrebări și cu ajutorul limbajului R, prin simulare. Din abordarea frecvenționistă am văzut că prin repetarea de \(N\) ori a unui experiment în condiții identice,
\[ \mathbb{P}(A|B)\approx \frac{N(A\cap B)}{N(B)} \]
unde \(N(A\cap B)\) este numărul de realizări (din \(N\)) a evenimentului \(A\cap B\) iar \(N(B)\) este numărul de realizări a evenimentului \(B\). Să considerăm \(N = 10^5\) și fie
N <- 10^5
copil1 <- sample(c("baiat", "fata"), N, replace = TRUE)
copil2 <- sample(c("baiat", "fata"), N, replace = TRUE)Aici copil1 este un vector de lungime N care reprezintă sexul primului copil și în mod similar copil2 reprezintă sexul celui de-al doilea copil.
Fie \(A\) evenimentul ca ambii copii să fie băieți și \(B\) evenimentul prin care cel mai tânăr este băiat.
nB <- sum(copil2 == "baiat")
nAB <- sum(copil1 == "baiat" & copil2 == "baiat")
p2 <- nAB/nBPrin urmare probabilitatea (simulată) ca familia să aibă cei doi copii băieți știind că cel tânăr este băiat este 0.4985668.
Considerând acum \(C\) evenimentul prin care familia are cel puțin un copil băiat, avem
nC <- sum(copil1 == "baiat" | copil2 == "baiat")
p1 <- nAB/nCde unde probabilitatea ca familia să aibă doar băieți știind că cel puțin unul este băiat este 0.3322913.
Să presupunem acum că vrem să rezolvăm următoarea problemă
Exercițiul 1.2 O familie are doi copii. Să se determine probabilitatea ca ambii copii să fie de sex feminin știind că cel puțin unul dintre ei este fata născută iarna. Vom presupune că avem șanse egale ca un copil să se fi născut în oricare din cele patru anotimpuri și că între sexul copilului și anotimp nu există nicio legătură.
Soluție. Pentru a răspunde la întrebare, vom defini spațiul stărilor prin
\[ \Omega = \{BP,BV,BT,BI, FP, FV, FT, FI\}^2 \]
unde \(B, F\) înseamnă băiat, respectiv fată iar \(P, V, T, I\) fac referire la anotimp (primăvară, vară, toamnă și respectiv iarnă), astfel un element \(\omega\in\Omega\) este de forma \(\omega=(BT, FI)\) și se citește primul copil este de sex masculin născut toamna iar cel de-al doilea copil este de sex feminin născut iarna. În ipotezele suplimentare din problemă, ne aflăm în contextul câmpului de probabilitate al lui Laplace, în care \(\mathbb{P}(\omega)=\frac{1}{64}\) pentru orice \(\omega\in\Omega\).
Problema ne cere să calculăm
\[ \mathbb{P}(\text{ambii să fie F } | \text{ cel puțin unul este FI} ) = ? \]
Cum
\[ \mathbb{P}(\text{ cel puțin unul este FI} ) = \frac{\#\{\text{cel puțin unul este FI}\}}{64} \]
iar \(|\{\text{cel puțin unul este FI}\}| = 15\), deoarece sunt \(8\) elemente de forma \((FI,\cdot)\) (i.e. care încep cu \(FI\)) și 7 elemente de forma \((\cdot,FI)\) (i.e. care se termină cu \(FI\)) găsim că
\[ \mathbb{P}(\text{ cel puțin unul este FI} ) = \frac{15}{64}. \]
În mod similar, probabilitatea
\[ \mathbb{P}(\text{ambii să fie F } \cap \text{ cel puțin unul este FI} ) = \frac{\#\{\text{ambii să fie F și cel puțin unul este FI}\}}{64} = \frac{7}{64} \]
deoarece sunt \(4\) elemente de forma \((FI, F\_ )\) și \(3\) de forma \((F\_, FI)\), fără să numărăm de două ori \((FI,FI)\).
Astfel, probabilitatea căutată este
\[ \mathbb{P}(\text{ambii să fie F } | \text{ cel puțin unul este FI} ) = \frac{7}{15}. \]
Ca și în cazul soluției din Exercițiul 1.1, putem găsi răspunsul la problemă și prin simulare. Mai precis, generând \(N = 10^5\) familii, avem
# generam N familii
N <- 10000 # numarul de repetari ale experimentului
copil1 <- sample(c("baiat-primavara", "baiat-vara", "baiat-toamna", "baiat-iarna",
"fata-primavara", "fata-vara", "fata-toamna", "fata-iarna"), N, replace = TRUE)
copil2 <- sample(c("baiat-primavara", "baiat-vara", "baiat-toamna", "baiat-iarna",
"fata-primavara", "fata-vara", "fata-toamna", "fata-iarna"), N, replace = TRUE)Dacă notăm cu \(B=\{\text{cel puțin unul este FI}\}\) evenimentul cu care condiționăm, atunci frecvența de apariție acestuia este
nB <- sum(copil1 == "fata-iarna" | copil2 == "fata-iarna")iar dacă \(A=\{\text{ambii să fie F}\}\), frecvența de apariție a evenimentului \(A\cap B\) este
nAB <- sum(copil1 == "fata-iarna" & (copil2 == "fata-vara" | copil2 == "fata-primavara" | copil2 == "fata-toamna" | copil2 == "fata-iarna")) +
sum(copil2 == "fata-iarna" & (copil1 == "fata-vara" | copil1 == "fata-primavara" | copil1 == "fata-toamna"))Probabilitatea căutată, \(\mathbb{P}(A|B)\), o putem estima prin
p_AB <- nAB/nBși obținem 0.4814657 (probabilitatea teoretică fiind \(\frac{7}{15}=\) 0.4666667).
Aruncatul cu zarul
Exercițiul 1.3 Să presupunem că aruncăm trei zaruri echilibrate (sau cu un zar de trei ori). Care este probabilitatea ca suma punctelor de pe fața superioară să fie egal cu \(8\)? Care este probabilitatea ca primul zar (sau la prima aruncare) să arate mai puțin de \(3\) puncte dat fiind că suma punctelor de pe fețele celor trei zaruri este mai mare sau egală cu \(8\)?
Soluție. Vom considera următoarea funcție în R:
prob_zar <- function(d = 3, s = 8, n = 1000){
# d - numarul de zaruri
# s - suma valorilor celor d zaruri
# n - nr de repetari ale experimentului
rez <- rep(0, n)
# repet experimentul de n ori
for (i in 1:n){
x <- sample(1:6, d, replace = TRUE)
rez[i] <- sum(x)
}
return(list(r1 = sum(rez == s)/n,
r2 = table(rez)/n))
}Obținem că pentru \(n = 10000\) de repetări ale experimentului, probabilitatea este:
z <- prob_zar(3, 8, 10000)
z$r1[1] 0.1016Următoarea figură ilustrează repartiția sumei celor trei zaruri:
barplot(height = as.numeric(z$r2),
names.arg = names(z$r2))
Pentru o scriere mai compactă putem folosi funcția replicate:
zard <- function(d){
return(sample(1:6, d, replace = TRUE))
}
prob_zar2 <- function(d = 3, s = 8, n = 1000){
rez <- replicate(n, sum(rolld(d)))
return(list(r1 = sum(rez == s)/n,
r2 = table(rez)/n))
}Pentru partea a doua a exercițiului avem de-a face cu o probabilitate condiționată. Dacă notăm cu \(A\) evenimentul prin care suma zarurilor este \(\geq 8\) și cu \(B\) evenimentul prin care primul zar are mai puțin de \(3\) puncte atunci căutăm \(\mathbb{P}(B|A)\). Următorul cod ne permite aproximarea acestei probabilități:
n <- 1000
cA <- 0 # suma zarurilor
cBA <- 0 # X1 < 3 si suma >= 8
for (i in 1:n){
x <- sample(1:6, 3, replace = TRUE)
if (sum(x) >= 8){
cA <- cA + 1
if (x[1] < 3){
cBA <- cBA + 1
}
}
}
# estimarea probabilitatii
cBA/cA[1] 0.2646007Monty Hall
Exercițiul 1.4 Sunteți participant într-un joc televizat în care gazda vă prezintă trei uși închise. Acesta vă spune că în spatele unei uși se află o mașină iar în spatele celorlalte două se află câte o capră. Jocul decurge în felul următor: trebuie să alegeți una dintre cele trei uși; gazda, care știe în spatele cărei uși se află mașina, deschide una dintre celelalte două uși, în spatele căreia se află o capră apoi vă întreabă dacă vreți să rămâneți la alegerea inițială sau vreți să alegeți cealaltă ușă rămasă închisă. Presupunând că vreți să câștigați o mașină, ce alegere preferați ?
Soluție. Considerăm evenimentele \(U_1\), \(U_2\) și \(U_3\) ca indicând ușa în spatele căreia se află mașina. Aceste evenimente au probabilitatea de \(1/3\) fiecare.
Evenimentul prin care jucătorul alege ușa cu numărul 1 este notat cu \(J_1\) și evenimentul prin care gazda deschide ușa cu numărul 3 este notat cu \(G_3\). Avem că \(\mathbb{P}(U_i|J_1) = 1/3\) iar \(\mathbb{P}(G_3|U_1,J_1) = 1/2\), \(\mathbb{P}(G_3|U_2,J_1) = 1\) și \(\mathbb{P}(G_3|U_3,J_1) = 0\).
Avem că probabilitatea ca jucătorul să câștige adoptând schimbarea ușilor este (în condițiile în care a ales inițial poarta 1 și gazda a deschis poarta 3)
\[ \begin{aligned} \mathbb{P}(U_2|G_3,J_1) &= \frac{\mathbb{P}(G_3|U_2,J_1)\mathbb{P}(U_2\cap J_1)}{\mathbb{P}(G_3\cap J_1)}\\ &=\frac{\mathbb{P}(G_3|U_2,J_1)\mathbb{P}(U_2\cap J_1)}{\mathbb{P}(G_3|U_1,J_1)\mathbb{P}(U_1\cap J_1)+\mathbb{P}(G_3|U_2,J_1)\mathbb{P}(U_2\cap J_1)+\mathbb{P}(G_3|U_3,J_1)\mathbb{P}(U_3\cap J_1)}\\ &= \frac{\mathbb{P}(G_3|U_2,J_1)}{\mathbb{P}(G_3|U_1,J_1)+\mathbb{P}(G_3|U_2,J_1)+\mathbb{P}(G_3|U_3,J_1)} = \frac{1}{1/2+1+0} = \frac{2}{3} \end{aligned} \]
Intuiție (să presupunem că ne aflăm în situația în care am ales ușa cu numărul 1):
| Ușa 1 | Ușa 2 | Ușa 3 | Nu schimbăm | Schimbăm |
|---|---|---|---|---|
| Mașină | Capră | Capră | Pierdut | Câștigat |
| Capră | Mașină | Capră | Câștigat | Pierdut |
| Capră | Capră | Mașină | Câștigat | Pierdut |
Vom încerca să simulăm jocul descris de problema noastră:
monty <- function(random = TRUE) {
doors <- 1:3
if (random){
# alege aleator unde se afla masina
cardoor <- sample(doors,1)
}else{
print("Scrieti numarul usii in spatele careia se afla masina:")
cardoor = scan(what = integer(), nlines = 1, quiet = TRUE)
}
# Monty ii spune jucatorului sa aleaga usa
print("Monty Hall spune ‘Alege o usa, orice usa!’")
# alegerea jucatorului (1,2 sau 3)
chosen <- scan(what = integer(), nlines = 1, quiet = TRUE)
# Monty intoarce o usa cu capra (nu poate fi usa aleasa
# de jucator si nici usa cu masina)
if (chosen != cardoor){
montydoor <- doors[-c(chosen, cardoor)]
}else{
montydoor <- sample(doors[-chosen],1)
}
# jucatorul schimba sau...
print(paste("Monty deschide usa ", montydoor, "!", sep=""))
print("Doresti sa schimbi usa (y/n)?")
reply <- scan(what = character(), nlines = 1, quiet = TRUE)
# ce incepe cu "y" este da
if (substr(reply,1,1) == "y"){
chosen <- doors[-c(chosen,montydoor)]
}
# Rezultatul jocului
if (chosen == cardoor){
print("Bravo! Ai castigat !")
}else{
print("Pacat! Ai pierdut!")
}
}Problema potrivirilor
Exercițiul 1.5 Să ne imaginăm că \(20\) de persoane merg la operă și că fiecare persoană poartă o pălărie. În momentul în care ajung la intrare își lasă pălăria la garderobă. Pe parcursul reprezentării artistice, persoana responsabilă cu garderoba se încurcă, pierde lista cu numerele locațiilor și returnează în mod aleator pălăriile persoanelor la plecare. Care este probabilitatea ca cel puțin o persoană să fi primit pălăria cu care a venit?
Soluție. Să presupunem că persoanele și pălăriilee sunt numerotate de la \(1\) la \(n = 20\) și inițial persoana \(i\) a venit cu pălăria \(i\). Fie \(\Omega = S_n\) mulțimea permutărilor cu \(n\) elemente, \(\mathcal{F} = \mathcal{P}(\Omega)\) și \(\mathbb{P}\) echiprobabilitatea pe \((\Omega, \mathcal{F})\). Să notăm cu \(E_i\) evenimentul prin care a \(i\)-a persoană primește pălăria cu numărul \(i\), adică primește pălăria cu care a venit. Evenimentul \(A\) prin care cel puțin o persoană a primit pălăria cu care a venit este
\[ A = E_1 \cup E_2 \cup\cdots\cup E_n. \]
Cum \(E_i\) reprezintă mulțimea permutărilor \(\sigma\in\S_n = \Omega\) pentru care \(\sigma(i) = i\) avem că
\[ \mathbb{P}(E_i) = \frac{(n-1)!}{n!} = \frac{1}{n} \]
deoarece \(|\Omega| = n!\) iar cele \(n-1\) valori diferite de \(i\) pot fi așezate în \((n-1)!\) moduri. În mod similar,
\[ \mathbb{P}(E_i\cap E_j) = \frac{(n-2)!}{n!} \]
și în general
\[ \mathbb{P}(E_{i_1}\cap E_{i_2}\cap\cdots\cap E_{i_k}) = \frac{(n-k)!}{n!}. \]
Formula lui Poincare permite calcularea probabilității dorite
\[ \mathbb{P}(A) = \mathbb{P}(E_1 \cup E_2 \cup\cdots\cup E_n) = \sum_{k = 1}^{n} (-1)^{k-1} \sum_{1\leq i_1<\cdots<i_k\leq n}\mathbb{P}(E_{i_1}\cap E_{i_2}\cap\cdots\cap E_{i_k}) \]
de unde
\[ \mathbb{P}(A) = \sum_{k = 1}^{n} (-1)^{k-1} \binom{n}{k}\frac{(n-k)!}{n!} = \sum_{k = 1}^{n} \frac{(-1)^{k-1}}{k!} \to 1 - \frac{1}{e} \approx 0.63212 \]
Să vedem că același rezultat îl obținem și prin simulare:
# repetam experimentul de un nr mare de ori
m <- 10000
# nr de persoane
n <- 20
n.potriviri <- rep(0, m)
persoane <- 1:n
for (i in 1:m){
palarii <- sample(1:n, n)
n.potriviri[i] <- sum(palarii == persoane)
}
# proportia persoanelor cu cel putin o potrivire
sum(n.potriviri > 0) / m[1] 0.6306Ruina jucătorului
Exercițiul 1.6 Un bărbat vrea să își cumpere un obiect (de exemplu o mașină sau o casă) care costă \(N\) unități monetare. Să presupunem că el are economisit un capital de \(0 < k < N\) unități monetare și încearcă să câștige restul jucând un joc de noroc cu managerul unei bănci. Jocul este următorul: bărbatul aruncă o monedă echilibrată în mod repetat iar dacă moneda pică cap (\(H\)) atunci managerul îi dă o unitate monetară, în caz contrar bărbatul plătește o unitate monetară băncii. Jocul continuă până când unul din două evenimente se realizează: sau bărbatul câștigă suma necesară și își cumpără obiectul dorit sau pierde banii și ajunge la faliment. Ne întrebăm care este probabilitatea să ajungă la faliment?
Soluție. Fie \(A\) evenimentul ca bărbatul să ajungă la ruină și \(B\) evenimentul ca la prima aruncare moneda a picat cap. Atunci din formula probabilității totale avem
\[ \mathbb{P}_k(A) = \mathbb{P}_k(A|B)\mathbb{P}(B)+\mathbb{P}_k(A|B^c)\mathbb{P}(B^c) \]
unde \(\mathbb{P}_k\) este probabilitatea calculată în funcție de valoarea \(k\) a capitalului inițial al jucătorului. Să observăm că \(\mathbb{P}_k(A|B)\) devine \(\mathbb{P}_{k+1}(A)\) deoarece dacă la prima aruncare avem cap atunci capitalul inițial a crescut la \(k+1\). În mod similar, dacă la prima aruncare am obținut coadă atunci \(\mathbb{P}_k(A|B^c) = \mathbb{P}_{k-1}(A)\). Notând cu \(p_k = \mathbb{P}_k(A|B)\) obținem următoarea ecuație
\[ p_k = \frac{1}{2}p_{k+1} + \frac{1}{2}p_{k-1}, \]
cu valorile inițiale \(p_0=1\) (dacă jucătorul a pornit cu un capital inițial nul atunci el este în faliment) și respectiv \(p_N=0\) (dacă jucătorul are din start suma necesară pentru a achiziționa obiectul dorit atunci nu mai are loc jocul).
O simulare a jocului pentru \(N = 50\) și \(k = 5\) este prezentată de următoarea funcție:
ruina <- function(N, k){
flag <- TRUE
joc <- 0
capital <- k
y <- capital
while(flag){
x <- 2*rbinom(1,1,0.5)-1
capital <- capital + x
y <- c(y, capital)
joc <- joc + 1
if (capital == 0 || capital == N){
flag = FALSE
}
}
return(y) # daca am 0 este ruina altfel este succes
}Putem ilustra grafic jocul după cum urmează:
N <- 50
k <- 5
set.seed(1234)
y <- ruina(N, k)
joc <- length(y) - 1 # nr de jocuri
plot(0:joc, y, type = "l",
main = "Ruina jucatorului",
xlab = "Numar de jocuri",
ylab = "Capitalul jucaturului",
bty = "n",
lty = 3, col = "grey50")
abline(h = c(0,N), col = "lightgrey", lty = 3)
points(0:joc, y, col = myblue,
pch = 16,
cex = 0.5)
points(0, k, col = myred, pch = 16)
text(0, k, labels = paste0("k = ", k), pos = 1, cex = 0.7)
Dacă definim \(b_k = p_k - p_{k-1}\) pentru \(k\geq 1\) atunci \(b_k = b_{k-1}\), \(\forall k\geq2\). Prin urmare \(b_k = b_1\) și \(p_k = b_1+p_{k-1} = kb_1+p_0\). Observând că \(b_1+\cdots+b_N=p_N-p_0=-1\) deducem că \(b_1=-\frac{1}{N}\) iar \(p_k=1-\frac{k}{N}\).
Dorim să repetăm experimentul de \(M = 1000\) de ori (pentru valorile inițiale \(N = 50\) și \(k = 5\)) și ne interesăm de câte ori jucătorul a ajuns la faliment.
N <- 50
k <- 5
M <- 1000
# Obs. - rezultatul functiei ruina trebuie modificat
joc <- replicate(M, ruina(N, k)) # repeta functia de M ori
proba_ruina <- sum(joc == 0)/M Am obținut că probabilitatea empirică de faliment este 0.898 iar cea teoretică este 0.9.
Aruncatul cu banul (II)
Exercițiul 1.7 O monedă are probabilitatea să pice cap \(p\) și să pice coadă \(q\) astfel ca \(p+q=1\). Moneda este aruncată succesiv și independent până când evenimentul \(A\), am obținut două capete unul după altul sau două cozi una după alta, s-a realizat. Determinați numărul mediu de aruncări necesare realizării evenimentului de interes, \(A\).
Soluție. Vom prezenta două soluții pentru acest exercițiu. În prima soluție vom începe prin a calcula funcția de masă pentru variabila care descrie experimentul și apoi vom calcula media. În cea de-a doua soluție vom calcula media cu ajutorul condiționării.
Fie \(X\) variabila aleatoare care reprezintă numărul de aruncări necesare pentru ca evenimentul \(A\) să se realizeze (numărul de aruncări necesare obținerii a \(2\) capete unul după altul sau a \(2\) cozi una după alta) și \(X_i\in\{H,T\}\) variabila aleatoare care descrie rezultatul obținut la cea de-a i-a aruncare. Evenimentul \(\{X=n\}\) poate fi exprimat ca reuniunea dintre \(A_n\) și \(B_n\), unde \(A_n\) reprezintă evenimentul să obținem două capete consecutive pentru prima oară la cea de-a \(n\)-a aruncare iar \(B_n\) este evenimentul să obținem două cozi consecutive la a \(n\)-a aruncare.
Se poate observa cu ușurință că un eveniment elementar din \(A_n\) are forma \(\omega = \underbrace{\cdots}_{n-2} HH\) și se poate determina în întregime datorită faptului că în primele \(n-2\) aruncări nu putem avea două realizări consecutive. În plus dacă \(n = 2k+1\) atunci
\[ \mathbb{P}(A_{2k+1}) = \mathbb{P}(X_1 = T, X_2 = H, \cdots, X_{2k-1}=T, X_{2k} = H, X_{2k+1} = H) \overset{indep.}{=} (pq)^kp \]
iar dacă \(n=2k+2\) atunci
\[ \mathbb{P}(A_{2k+2}) = \mathbb{P}(X_1 = H, X_2 = T, \cdots, X_{2k}=T, X_{2k+1} = H, X_{2k+2} = H) \overset{indep.}{=} (pq)^kp^2 \]
În mod similar se poate calcula probabilitatea evenimentului \(B_n\) pentru \(n=2k+1\)
\[ \mathbb{P}(B_{2k+1}) = \mathbb{P}(X_1 = H, X_2 = T, \cdots, X_{2k-1} = H, X_{2k} = T, X_{2k+1} = T) \overset{indep.}{=} (pq)^kq \]
și respectiv \(n=2k+2\)
\[ \mathbb{P}(B_{2k+2}) = \mathbb{P}(X_1 = T, X_2 = H, \cdots, X_{2k} = H, X_{2k+1} = T, X_{2k+2} = T) \overset{indep.}{=} (pq)^kq^2 \]
Cum \(\mathbb{P}(X=n) = \mathbb{P}(A_n) + \mathbb{P}(B_n)\) deducem că
\[ \mathbb{P}(X=n) =\left\{\begin{array}{ll} (pq)^{\frac{n-1}{2}}, & \text{$n$ impar}\\ (pq)^{\frac{n-2}{2}}(p^2+q^2), & \text{$n$ par} \end{array}\right. \]
Următoarea funcție simulează evenimentul din enunțul problemei:
flip_coins <- function(p){
flip <- sample(c("H", "T"), 1, prob = c(p, 1-p))
nflips <- 1
flag <- TRUE
while(flag){
x <- sample(c("H", "T"), 1, prob = c(p, 1-p))
nflips <- nflips + 1
if (flip == x){
flag <- FALSE
}
flip <- x
}
return(nflips)
}
flip_coins(0.2)[1] 3Înainte de a calcula media ne propunem să repetăm de \(N = 10000\) de ori experimentul (pentru \(p=0.3\)) și să comparăm rezultatul empiric cu cel teoretic.
N <- 10000
p <- 0.3
rez_flips <- replicate(N, flip_coins(p))Rezultatele sunt incluse în tabelul de mai jos:
| n | Empiric | Teoretic |
|---|---|---|
| 2 | 0.5760 | 0.58000 |
| 3 | 0.2110 | 0.21000 |
| 4 | 0.1245 | 0.12180 |
| 5 | 0.0434 | 0.04410 |
| 6 | 0.0255 | 0.02558 |
| 7 | 0.0100 | 0.00926 |
| 8 | 0.0059 | 0.00537 |
| 9 | 0.0016 | 0.00194 |
| 10 | 0.0014 | 0.00113 |
| 11 | 0.0004 | 0.00041 |
| 12 | 0.0002 | 0.00024 |
| 13 | 0.0001 | 0.00009 |
Pentru a calcula media avem
\[ \mathbb{E}[X] = \sum_{n=2}^{\infty}n\mathbb{P}(X=n) \]
iar dacă seriile \(\sum_{k=1}^{\infty}(2k+1)\mathbb{P}(X=2k+1)\) și \(\sum_{k=0}^{\infty}(2k+2)\mathbb{P}(X=2k+2)\) sunt convergente atunci putem scrie
\[ \mathbb{E}[X] = \sum_{n=2}^{\infty}n\mathbb{P}(X=n) = \sum_{k=1}^{\infty}(2k+1)\mathbb{P}(X=2k+1) + \sum_{k=0}^{\infty}(2k+2)\mathbb{P}(X=2k+2). \]
Se poate arăta cu ușurință că
\[ \sum_{k=1}^{\infty}(2k+1)\mathbb{P}(X=2k+1) = 2\sum_{k=1}^{\infty}(k+1)(pq)^k - \sum_{k=1}^{\infty}(pq)^k = \frac{pq(3-pq)}{(1-pq)^2} \]
și
\[ \sum_{k=0}^{\infty}(2k+2)\mathbb{P}(X=2k+2) = 2(p^2+q^2)\sum_{k=0}^{\infty}(k+1)(pq)^k = \frac{2(p^2+q^2)}{(1-pq)^2} \]
de unde concluzionăm că \(\mathbb{E}[X] = \frac{pq(3-pq)}{(1-pq)^2} + \frac{2(p^2+q^2)}{(1-pq)^2} = \frac{2+pq}{1-pq}\).
O a doua soluție pentru determinarea mediei este bazată pe condiționare. Fie \(H_k\) și respectiv \(T_k\) evenimentele prin care capul respectiv coada a apărut la cea de-a \(k\)-a aruncare și \(\mathbb{P}(H_k) = p\) iar \(\mathbb{P}(T_k) = q\). Din formula probabilității totale, avem că
\[ \mathbb{E}[X] = \mathbb{E}[X|H_1]\mathbb{P}(H_1) + \mathbb{E}[X|T_1]\mathbb{P}(T_1) = p\mathbb{E}[X|H_1] + q\mathbb{E}[X|T_1]. \]
Mai mult,
\[ \mathbb{E}[X|H_1] = p\mathbb{E}[X|H_1\cap H_2] + q\mathbb{E}[X|H_1\cap T_2] \]
și cum \(\mathbb{E}[X|H_1\cap H_2] = 2\) (evenimentul \(A\) s-a realizat la a doua aruncare) iar \(\mathbb{E}[X|H_1\cap T_2] = 1 + \mathbb{E}[X|T_1]\) (dacă jocul nu este gata doar ultima aruncare este importantă) obținem că
\[ \mathbb{E}[X|H_1] = 2p + q(1+\mathbb{E}[X|T_1]). \]
În mod similar găsim \(\mathbb{E}[X|T_1] = 2q + p(1+\mathbb{E}[X|H_1])\). Rezolvând sistemul de două ecuații cu două necunoscute obținem soluțiile \(\mathbb{E}[X|H_1] = \frac{2+q^2}{1-pq}\) și \(\mathbb{E}[X|T_1] = \frac{2+p^2}{1-pq}\) de unde media este \(\mathbb{E}[X] = \frac{2+pq}{1-pq}\).
Putem verifica numeric dacă formula pe care am gasit-o este corectă. Pentru aceasta vom repeta experimentul (\(p = 0.3\)) de \(N = 10000\) de ori.
Obținem că media empirică este 2.807 iar cea teoretică este 2.797.
Principiul reflexiei și problema scrutinului
Să presupunem că avem \(n\) elemente \(\epsilon_1, \ldots, \epsilon_n\in\{-1,+1\}\) astfel încât \(p\) dintre ele iau valori de \(+1\) și \(q\) dintre ele iau valori de \(-1\) (\(n = p + q\)). Sumele parțiale \(s_k = \epsilon_1+\epsilon_2+\cdots+\epsilon_k\) reprezintă diferența dintre numărul de elemente de \(+1\) și de elemente de \(-1\) în primele \(k\) elemente. Observăm că
\[ s_k-s_{k-1}=\epsilon_k=\pm 1,\quad s_0 = 0,\quad s_n = p-q. \]
Din punct de vedere geometric, \(n\)-uplul \((\epsilon_1,\ldots,\epsilon_n)\) poate fi reprezentat cu ajutorul unei linii poligonale ce pleacă din origine și ajunge în punctul de coordonate \((n,s_n)\), panta fiecărui segment de lungime \(1\) este dată de \(\epsilon_k\) (\(+1\) NE și \(-1\) SE). Linia poligonală \((s_1,\ldots,s_n)\) are al \(k\)-lea punct de abscisă \(k\) și ordonată \(s_k\).
Dacă \(n\) și \(x\) sunt numere naturale nenule astfel ca \(n = p+q\) și \(x = p-q\), \(p, q\in\mathbb{N}\) atunci numărul de linii poligonale (în sensul de mai sus) din origine către punctul \((n,x)\) este dat de
\[ N_{n,x} = \binom{p+q}{p}. \]
Exercițiul 1.8 (Problema scrutinului) Să presupunem că în turul doi al alegerilor prezidențiale din 2019 participă doi candidați, \(P\) și \(Q\). Dacă voturile sunt date de manieră independentă cu probabilitatea de \(\frac{1}{2}\) pentru fiecare candidat și candidatul \(P\) primește \(p\) voturi iar candidatul \(Q\) primește \(q\) voturi astfel ca \(P\) să câștige (\(p>q\)), atunci să se calculeze probabilitatea ca pe tot parcursul numărării voturilor candidatul \(P\) să fi avut mai multe voturi decât candidatul \(Q\).
Soluție. Problema scrutinului (propusă și rezolvată de Whitworth în 18781) poate fi interpretată geometric prin intermediul unei linii poligonale de lungime \(p+q\) în care pantele fiecărui segment reprezintă opțiunea votului pentru unul din cei doi candidați: \(\epsilon_k = +1\) dacă al \(k\)-lea vot a fost pentru candidatul \(P\) și \(\epsilon_k = -1\) dacă al \(k\)-lea vot a fost pentru candidatul \(Q\). În mod similar, dată fiind o linie poligonală care pleacă din origine și ajunge în punctul \((p+q, p-q)\), aceasta poate reprezenta parcursul unui scrutin în care cei doi candidați vor avea la final \(p\) și respectiv \(q\) voturi.
Conform notațiilor precedente, putem observa că \(s_k\) reprezintă numărul voturilor cu care candidatul \(P\) conduce (sau este în urmă) imediat după cel de-al \(k\)-lea vot. Astfel cerința problemei se poate traduce în modul următor: candidatul \(P\) conduce pe tot parcursul procesului de votare dacă și numai dacă \(s_1>0, s_2>0,\ldots, s_n>0\) (linia poligonală se află tot timpul deasupra axei absciselor).
Soluția problemei se bazează pe următorul rezultat, numit principiul reflexiei. Acesta spune că numărul de traiectorii (linii poligonale în sensul descris mai sus) de la punctul de coordonate \((x,y)\) la punctul de coordonate \((x',y')\), cu \(x'>x\geq0\), \(y>0\) și \(y'>0\), care intersectează axa absciselor este egal cu numărul de traiectorii de la punctul de coordonate \((x,-y)\) la punctul de coordonate \((x',y')\).
Pentru a verifica această afirmație, fie \((s_x = y, s_{x+1},\ldots,s_{x'} = y')\) o linie poligonală (\(s_k-s_{k-1}=\pm 1\), \(k\in\{x,x+1,\ldots, x'\}\)) de la punctul \((x,y)\) la punctul \((x',y')\) care să intersecteze axa absciselor și fie \(t\) abscisa primului punct de intersecție (alegem \(t\) astfel încât \(s_x>0,\ldots,s_{t-1}>0\) și \(s_t = 0\) - a se vedea figura). Prin urmare \((-s_x = -y, -s_{x+1},\ldots,-s_{t-1},0,s_{t+1},\ldots,s_{x'} = y')\) este o linie poligonală de la punctul \((x,-y)\) la punctul \((x',y')\) care intersectează axa absciselor pentru prima oară în punctul \((t,0)\). Cum secțiunea \((s_x, s_{x+1},\ldots,s_{t-1},0)\) este reflexia secțiunii \((-s_x, -s_{x+1},\ldots,-s_{t-1},0)\), deducem că putem construi o bijecție de la mulțimea liniilor poligonale dintre \((x,y)\) și \((x',y')\) la mulțimea liniilor poligonale dintre \((x,-y)\) și \((x',y')\).
Vom arăta că numărul de linii poligonale care pleacă din origine și ajung în punctul de coordonate \((p+q,p-q)\) verificând \(s_1>0, s_2>0,\ldots, s_n>0\) este dat de formula:
\[ N = \frac{p-q}{p+q}\binom{p+q}{p}. \]
Pentru a verifica această formulă să observăm că numărul de linii poligonale de la punctul \((0,0)\) la punctul \((n,x)\) (\(n = p+q\) și \(x = p-q\)) care verifică \(s_1>0, s_2>0,\ldots, s_n>0\) coincide cu numărul de traiectorii de la \((1,1)\) la \((n,x)\). Acest număr reprezintă diferența dintre numărul de traiectorii din origine la punctul \((n-1,x-1)\) (translatăm toate punctele spre stânga și în jos cu o unitate), \(N_{n-1,x-1}\), mai puțin acele traiectorii care intersectează axa absciselor. Folosind principiul reflexiei, numărul de traiectorii care pornesc din \((1,1)\) și ajung în \((n,x)\) intersectând axa absciselor este egal cu numărul de traiectorii care pleacă din \((1,-1)\) și ajung în \((n,x)\), care este egal cu \(N_{n-1,x+1}\) (translatăm spre stânga și în sus cu o unitate). Prin urmare avem
\[ N = N_{n-1,x-1} - N_{n-1,x+1} = \binom{p+q-1}{p-1} - \binom{p+q-1}{p} = \frac{p-q}{p+q}\binom{p+q}{p}. \]
Probabilitatea căutată devine astfel \(\frac{N_{n-1,x-1} - N_{n-1,x+1}}{N_{n,x}} = \frac{p-q}{p+q}\).
Pentru a verifica numeric rezultatul obținut să considerăm următoarea funcție care simulează scrutin-ul pentru cei doi candidați:
scrutin <- function(sd){
set.seed(sd)
# generam voturile
vectscrutin = sample(c(rep("P",p),rep("Q",q)))
# transformam in +-1
vectscrutin = 2*(vectscrutin == "P") - 1
# calculam s_k
s <- cumsum(vectscrutin)
# returnam daca avem sau nu indeplinita conditia
sum(s>0) == p+q
}Considerând \(p = 50\) și \(q = 35\) și repetând experimentul de \(M = 100000\) de ori obținem că probabilitatea empirică este 0.1767 iar cea teoretică este 0.1765.
Legea arcsinus-ului
Exercițiul 1.9 (Repartiția arcsinus) Să presupunem că aruncăm cu banul de \(100\) de ori și că înregistrăm la fiecare aruncare numărul de capete și numărul de cozi. Definim ultima egalitate ca fiind ultima aruncare în care numărul de capete este egal cu numărul de cozi (deci ia valori de la \(0\) la \(100\)). Ne propunem să scriem un program care să determine locația ultimei egalități în jocul de noroc descris.
Soluție. Vom începe prin a transpune problema în limbaj matematic2. Fie \((X_n)_n\) un șir de variabile aleatoare independente ce iau valori în mulțimea \(\{+1,-1\}\) cu probabilitatea \(p = \frac{1}{2}\) (aruncăm cu banul și ne deplasăm la dreapta sau la stânga în funcție de rezultatul aruncării) și fie \(S_n = X_1 + \cdots + X_n\). Putem interpreta șirul \(S_1,S_2,\ldots,S_n\) geometric, ca și în cazul problemei scrutinului, cu ajutorul unei linii poligonale cu segmentele \((k-1, S_{k-1})\to(k, S_k)\). Fie \(L_{2n}\) variabila aleatoare care ne dă timpul ultimei întoarceri în origine a unui mers la întâmplare de lungime \(2n\) (de ce avem \(2n\)?), cu alte cuvinte
\[ L_{2n} = \sup\{m\leq 2n\,|\,S_m = 0\}. \]
Variabila aleatoare \(L_{2n}\) este cea care descrie ultima egalitate din enunțul problemei și ne propunem să găsim repartiția acesteia. Observăm că
\[\begin{align*} \mathbb{P}(L_{2n} = 2k) &= \mathbb{P}(S_{2k} = 0, S_{2k+1}\neq 0, \ldots, S_{2n}\neq 0)\\ &= \mathbb{P}(S_{2k} = 0)\mathbb{P}(S_{2k+1}\neq 0, \ldots, S_{2n}\neq 0|S_{2k} = 0). \end{align*}\]
O realizare a variabilei \(L_{100}\), în contextul problemei atunci când \(n = 50\) este ilustrată în figura următoare.
Am văzut, în problema scrutinnului, că \(\mathbb{P}(S_n=m)\) este \(\frac{N_{n,m}}{2^n}\) (numărul de linii poligonale din origine până în punctul \((n,m)\) supra numărul total de linii poligonale cu \(n\) puncte) de unde
\[ \mathbb{P}(S_n=m) = \binom{n}{\frac{n+m}{2}}2^{-n}, \]
unde folosim convenția că \(\binom{n}{\frac{n+m}{2}} = 0\) dacă \(\frac{n+m}{2}\) nu este un întreg între \(0\) și \(n\). Astfel deducem că \(\mathbb{P}(S_{2n} = 0) = \binom{2n}{n}2^{-2n}\).
Vom arăta că pentru un mers la întâmplare care începe din origine are loc relația
\[ \mathbb{P}(S_1\neq 0, S_2\neq 0, \ldots, S_{2n}\neq 0) = \mathbb{P}(S_{2n}=0) \]
și cum \(X_i\) sunt independente (fiecare aruncare se face de manieră independentă) atunci probabilitatea condiționată \(\mathbb{P}(S_{2k+1}\neq 0, \ldots, S_{2n}\neq 0|S_{2k} = 0)\) este egală cu probabilitatea ca un mers la întâmplare de lungime \(2n-2k\), care pornește din origine, să nu mai viziteze niciodată originea în cei \(2n-2k\) pași sau
\[ \mathbb{P}(S_{2k+1}\neq 0, \ldots, S_{2n}\neq 0|S_{2k} = 0) = \mathbb{P}(S_{2n-2k}=0). \]
Prin urmare avem relația \(\mathbb{P}(L_{2n} = 2k) = \mathbb{P}(S_{2k}=0)\mathbb{P}(S_{2n-2k}=0)\). Pentru ca demonstrația să fie completă rămâne să verificăm că
\[ \mathbb{P}(S_1\neq 0, S_2\neq 0, \ldots, S_{2n}\neq 0) = \mathbb{P}(S_{2n}=0). \]
Cum originea nu mai este vizitată în cei \(2n\) pași atunci sau \(S_j>0\) sau \(S_j<0\) pentru toți \(j\in\{1,2,\ldots,n\}\). Din simetrie, aceste alternative au aceeași probabilitate deci
\[ \mathbb{P}(S_1\neq 0, S_2\neq 0, \ldots, S_{2n}\neq 0) = 2 \mathbb{P}(S_1> 0, S_2> 0, \ldots, S_{2n}> 0) \]
și aplicând problema scrutinului (linia poligonală se afla deasupra axei absciselor) avem
\[\begin{align*} \mathbb{P}(S_1> 0, S_2> 0, \ldots, S_{2n}> 0) &= \sum_{k = 1}^{n}\mathbb{P}(S_1> 0, S_2> 0, \ldots, S_{2n}> 0, S_{2n} = 2k)\\ &= \sum_{k = 1}^{n}\frac{1}{2^{2n}}(N_{2n-1, 2k-1} - N_{2n-1, 2k+1})\\ &= \frac{1}{2^{2n}}N_{2n-1, 1} = 2^{-2n}\binom{2n-1}{n} = 2^{-2n-1}\binom{2n}{n} \\ &= \frac{1}{2}\mathbb{P}(S_{2n}=0), \end{align*}\]
de unde deducem concluzia.
Se poate verifica cu ușurință că \(\mathbb{P}(S_{2m}=0) \approx \frac{1}{\sqrt{\pi m}}\)3 pentru \(m\) suficient de mare, astfel
\[ \mathbb{P}(L_{2n} = 2k) = \mathbb{P}(S_{2k}=0)\mathbb{P}(S_{2n-2k}=0) \approx \frac{1}{\sqrt{\pi n}}\times\frac{1}{\sqrt{\pi(n-k)}} = \frac{1}{\pi\sqrt{k(n-k)}}. \]
Observăm că pentru \(\frac{k}{n}\to x\) avem \(n\mathbb{P}(L_{2n} = 2k)\to \frac{1}{\pi\sqrt{x(1-x)}}\). Dacă \(0<a<b<1\) și considerăm \(2na_n\) cel mai mic întreg par mai mare decât \(2na\) și respectiv \(2nb_n\) cel mai mare întreg par mai mic decât \(2nb\) atunci
\[\begin{align*} \mathbb{P}\left(a < \frac{L_{2n}}{2n} < b\right) &= \sum_{k = na_n}{nb_n}\mathbb{P}(L_{2n} = 2k) \\ &\approx \sum_{k = na_n}{nb_n}\frac{1}{\pi\sqrt{k(n-k)}} \approx \frac{1}{\pi}\int_{na}^{nb}\frac{1}{\sqrt{x(n-x)}}\,dx \\ &\overset{y=nx}{=} \frac{1}{\pi}\int_{a}^{b}\frac{1}{\sqrt{y(1-y)}}\,dy. \end{align*}\]
În concluzie am obținut că
\[ n\mathbb{P}(L_{2n} = 2k) \to \frac{1}{\pi\sqrt{x(1-x)}} \]
și
\[ \mathbb{P}\left(a < \frac{L_{2n}}{2n} < b\right)\approx \frac{1}{\pi}\int_{a}^{b}\frac{1}{\sqrt{y(1-y)}}\,dy = \frac{2}{\pi}(\arcsin{\sqrt{b}} - \arcsin{\sqrt{a}}). \]
Următoarea histogramă, pentru care am considerat că experimentul constă din \(n=1000\) de aruncări cu banul și pe care l-am repetat de \(M = 10000\) de ori, ilustrează grafic legea arcsinusului pentru problema noastră:
Note de subsol
Soluția a fost publicată în cartea Choice and chance în 1886. În 1887 Joseph Bertrand a propus o formă mai generală a problemei care a fost rezolvată de către Désiré André cu ajutorul principiului reflexiei.↩︎
Această problemă este preluată din cartea lui W. Feller Introduction to probability and its applications, vol I, 1968, capitolul III.↩︎
Din formula lui Stirling avem că \(n!\approx \sqrt{2\pi}n^{n+\frac{1}{2}}e^{-n}\) prin urmare \(\binom{2n}{n}2^{-2n}\approx\frac{\sqrt{2\pi}(2n)^{2n+\frac{1}{2}}e^{-2n}}{\left(\sqrt{2\pi}n^{n+\frac{1}{2}}e^{-n}\right)^2}2^{-2n} = \frac{1}{\sqrt{\pi n}}\)↩︎