Probabilidade Básica em R

Conceitos Básicos

Um modelo probabilístico é uma descrição qualitativa de uma situação, fenômeno ou experimento cujo resultado é incerto. Para definir tal modelo seguimos os passos:

1. Descrevemos os possívies resultados do experimento (especificamos o espaço amostral)
2. Especificamos a lei probabilística, que atribui probabilidades aos resultados ou conjunto de resultados.

A probabilidade Pr(A) de um evento A é a proporção (aproximada) de vezes que o evento ocorre quando repetimos o experimento várias vezes, de maneira independente, sob as mesmas condições. Quanto maior o número de experimentos, melhor é esta aproximação.

Um evento é um resultado medido ou observado de um fenômeno, que aconteceu por chance.

Nas seções seguintes ilustraremos a ideia de probabilidade fazendo contagens de possibilidades e também por meio de experimentos computacionais usando o método de Monte Carlo.

Cartas de um baralho

Nesta seção utilizaremos mostraremos como utilizar o R para resolver problemas envolvendo objetos como cartas de um baralho. Embora os problemas listados abaixo possam resolvidos facilmente com uma calculadora, nosso objetivo é ilustrar o uso do R, com comandos que serão imprescidíveis em problemas bem mais complicados que abordaremos mais adiante. Para definir um baralho de cartas no R podemos proceder do seguinte modo, usando os comandos paste e expand.grid:

library(tidyverse)
library(gtools)

suits <- c("Diamonds", "Clubs", "Hearts", "Spades")
numbers <- c("Ace", "Deuce", "Three", "Four", "Five", "Six", "Seven", "Eight", "Nine", "Ten", "Jack", "Queen", "King")
deck <- expand.grid(number=numbers, suit=suits)
deck2 <- paste(deck$number,deck$suit)
deck2

 [1] "Ace Diamonds"   "Deuce Diamonds"
 [3] "Three Diamonds" "Four Diamonds" 
 [5] "Five Diamonds"  "Six Diamonds"  
 [7] "Seven Diamonds" "Eight Diamonds"
 [9] "Nine Diamonds"  "Ten Diamonds"  
[11] "Jack Diamonds"  "Queen Diamonds"
[13] "King Diamonds"  "Ace Clubs"     
[15] "Deuce Clubs"    "Three Clubs"   
[17] "Four Clubs"     "Five Clubs"    
[19] "Six Clubs"      "Seven Clubs"   
[21] "Eight Clubs"    "Nine Clubs"    
[23] "Ten Clubs"      "Jack Clubs"    
[25] "Queen Clubs"    "King Clubs"    
[27] "Ace Hearts"     "Deuce Hearts"  
[29] "Three Hearts"   "Four Hearts"   
[31] "Five Hearts"    "Six Hearts"    
[33] "Seven Hearts"   "Eight Hearts"  
[35] "Nine Hearts"    "Ten Hearts"    
[37] "Jack Hearts"    "Queen Hearts"  
[39] "King Hearts"    "Ace Spades"    
[41] "Deuce Spades"   "Three Spades"  
[43] "Four Spades"    "Five Spades"   
[45] "Six Spades"     "Seven Spades"  
[47] "Eight Spades"   "Nine Spades"   
[49] "Ten Spades"     "Jack Spades"   
[51] "Queen Spades"   "King Spades"   

Por exemplo, calculemos qual é a probabilidade de obtermos King ao retirarmos a primeira carta do baralho. Como um baralho tem 52 cartas e temos 4 Kings, um de cada naipe, a probabilidade de retirarmos uma carta com King é:

4/52

[1] 0.07692308

Podemos reobter este resultado a partir das cartas do baralho. Inicialmente determinamos as cartas com King:

kings <- paste("King", suits)
kings

[1] "King Diamonds" "King Clubs"    "King Hearts"  
[4] "King Spades"  

Examinamos, entre as cartas do baralho, quais são que tem King:

deck2 %in% kings

 [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [9] FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE
[17] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[25] FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE
[33] FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE
[41] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[49] FALSE FALSE FALSE  TRUE

O comando mean nos dá a proporção do número de resultados TRUE para o número de resultados FALSE:

mean( deck2 %in% kings)

[1] 0.07692308

que é aproximadamente 1/13, como esperado.

Calculemos agora qual é a probabilidade condicional de que a segunda carta seja King dado que a primeira carta for King. Como temos agora somente 3 cartas deste tipo e o baralho agora tem somente 51 cartas, esperamos que a probabilidade seja 3/51. Confirmemos este resultado em R.

Antes disso estudemos o comando permutations() que computa, para qualquer lista de tamanho n, todos diferentes modos que podemos selecionar r itens. Vejamos um exemplo: de quantos modos podemos selecionar 2 elementos de um grupo de 5 elementos?

permutations(5,2)

      [,1] [,2]
 [1,]    1    2
 [2,]    1    3
 [3,]    1    4
 [4,]    1    5
 [5,]    2    1
 [6,]    2    3
 [7,]    2    4
 [8,]    2    5
 [9,]    3    1
[10,]    3    2
[11,]    3    4
[12,]    3    5
[13,]    4    1
[14,]    4    2
[15,]    4    3
[16,]    4    5
[17,]    5    1
[18,]    5    2
[19,]    5    3
[20,]    5    4

Notemos que a ordem importa. Além disso números repetidos não ocorrem, pois uma vez que selecionamos um número, ele não aparece novamente. Opcionalmente o comando permutations pode ter um vetor como terceiro argumento. Por exemplo, façamos uma permutação (ou arranjo) de 4 números, 3,4,5,6 em grupos de 3 elementos. De outro modo, se temos uma caixa com 4 bolas marcadas com estes números, quais são os possíveis resultados ao selecionar 3 bolas aleatóriamente?

m <- permutations(4,3,c(3,4,5,6,7))
m

      [,1] [,2] [,3]
 [1,]    3    4    5
 [2,]    3    4    6
 [3,]    3    5    4
 [4,]    3    5    6
 [5,]    3    6    4
 [6,]    3    6    5
 [7,]    4    3    5
 [8,]    4    3    6
 [9,]    4    5    3
[10,]    4    5    6
[11,]    4    6    3
[12,]    4    6    5
[13,]    5    3    4
[14,]    5    3    6
[15,]    5    4    3
[16,]    5    4    6
[17,]    5    6    3
[18,]    5    6    4
[19,]    6    3    4
[20,]    6    3    5
[21,]    6    4    3
[22,]    6    4    5
[23,]    6    5    3
[24,]    6    5    4

Como vemos da matriz de saída do comando acima, o número de modos é 24, que pode ser obtido diretamente contando o número de linhas de m:

n <- nrow(m)
n

[1] 24

Façamos uma simulação de 6 escolhas aleatórias de 3 bolas

indice <- sample(n,6)
indice

[1] 13 23 20 15  2 16

Ou seja, as bolas selecionadas têm os números dados por:

m[indice,]

     [,1] [,2] [,3]
[1,]    5    3    4
[2,]    6    5    3
[3,]    6    3    5
[4,]    5    4    3
[5,]    3    4    6
[6,]    5    4    6

Exemplo:

Simule a seleção aleatória de de 5 números (0 a 9) de 7 dígitos (não mostre as possíveis seleções na tela).
Solução

numeros <- permutations(10,7,c(0:9))
n <- nrow(numeros)
n

[1] 604800

Selecionemos 5 números de 7 dígitos:

indice <- sample(n,5)
indice

[1] 149839  94593 349073 109609 408838

numeros[indice,]

     [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,]    2    5    3    4    7    9    6
[2,]    1    6    0    7    3    8    2
[3,]    5    7    9    3    8    4    0
[4,]    1    8    3    5    4    6    0
[5,]    6    7    8    4    9    5    1

Voltemos agora ao problema de calcular a probabilidade de que a segunda carta seja King dado que a primeira carta for King. Computemos todos os possíveis modos de escolher duas cartas, quando a ordem é relevante.

hands <- permutations(52,2,v=deck2)
nrow(hands)

[1] 2652

Ou seja, temos 2652 possibilidades. A primeira carta será selecionada da primeira coluna da matriz hands e similarmente com a segunda:

primeira_carta <- hands[,1]
segunda_carta <- hands[,2]

Vejamos em quantos casos a primeira carta é King:

k <- sum(primeira_carta %in% kings)
k

[1] 204

Ou seja, 204 casos. Para determinar a probabilidade condicional devemos agora determinar qual é a fração destes 204 casos em que temos King na segunda carta:

sum(primeira_carta %in% kings & segunda_carta %in% kings)/sum(primeira_carta %in% kings)

[1] 0.05882353

que é igual a 3/51, como esperado.

Notemos ainda que o mesmo resultado poderia ser obtido utilizando valores esperados:

mean(primeira_carta %in% kings & segunda_carta %in% kings)/mean(primeira_carta %in% kings)

[1] 0.05882353

Estas expressões acima são formas diferentes da relação fundamental \[ Pr(B|A)=\frac{Pr(A \cap B)}{Pr(A)} \] Suponhamos agora que a ordem não seja importante. No jogo blackjack se obtivermos um A e uma carta de face ou 10, teremos um 21 natural, a a vitória é instantânea. Queremos calcular a probabilidade disso acontecer. Notemos que agora devemos enumerar combinações e não permutações, pois a ordem não importa. Um A e um King ou um King e A são a mesma coisa e não queremos contar este caso duas vezes. Notemos a diferença entre permutaçôes, onde a ordem importa,

permutations(3,2)

     [,1] [,2]
[1,]    1    2
[2,]    1    3
[3,]    2    1
[4,]    2    3
[5,]    3    1
[6,]    3    2

e combinações, onde a ordem não importa:

combinations(3,2)

     [,1] [,2]
[1,]    1    2
[2,]    1    3
[3,]    2    3

Para calcular a probabilidade de um 21 natural, inicialmente definimos um vetor que contenha todos os A:

aces <- paste("Ace", suits)
aces

[1] "Ace Diamonds" "Ace Clubs"    "Ace Hearts"  
[4] "Ace Spades"  

Definimos agora um vetor que inclui todas as cartas de face e o 10:

facecard <- c("King","Queen", "Jack", "Ten")
facecard <- expand.grid(number=facecard, suit=suits)
facecard

Formamos o vetor explicitamente:

facecard <- paste(facecard$number,facecard$suit)
facecard

 [1] "King Diamonds"  "Queen Diamonds"
 [3] "Jack Diamonds"  "Ten Diamonds"  
 [5] "King Clubs"     "Queen Clubs"   
 [7] "Jack Clubs"     "Ten Clubs"     
 [9] "King Hearts"    "Queen Hearts"  
[11] "Jack Hearts"    "Ten Hearts"    
[13] "King Spades"    "Queen Spades"  
[15] "Jack Spades"    "Ten Spades"    

Calculamos todas as combinações de 2 cartas que podem ser obtidas a partir das 52 cartas:

hands <- combinations(52,2,v = deck2)
nrow(hands)

[1] 1326

Ou seja, temos 1326 pares de cartas.

summary(hands)

              V1                    V2      
 Ace Clubs     :  51   Three Spades  :  51  
 Ace Diamonds  :  50   Three Hearts  :  50  
 Ace Hearts    :  49   Three Diamonds:  49  
 Ace Spades    :  48   Three Clubs   :  48  
 Deuce Clubs   :  47   Ten Spades    :  47  
 Deuce Diamonds:  46   Ten Hearts    :  46  
 (Other)       :1035   (Other)       :1035  

Temos uma matriz com dois índices: o primeiro (1 a 1326) se refere ao ao par o segundo (1 a 2) à carta do par. Vejamos algumas cartas:

hands[1326,1]

[1] "Three Hearts"

hands[1326,2]

[1] "Three Spades"

hands[1325,1]

[1] "Three Diamonds"

hands[1,1]

[1] "Ace Clubs"

hands[1,2]

[1] "Ace Diamonds"

Por exemplo, o primeiro par é

hands[1,]

[1] "Ace Clubs"    "Ace Diamonds"

e o último par é:

hands[1326,]

[1] "Three Hearts" "Three Spades"

Calculemos agora as probabilidades:

mean(hands[,1] %in% aces & hands[,2] %in% facecard)

[1] 0.04826546

Notemos que embora tenhamos suposto acima que o A é a primeira carta, isso não é importante, pois nas combinações a ordem é irrelevante. Outro modo de calcular esta probabilidade é considerar ambas as possibilidades:

mean((hands[,1] %in% aces & hands[,2] %in% facecard) | (hands[,2] %in% aces & hands[,1] %in% facecard))

[1] 0.04826546

que dá o mesmo resultado.

Simulações de Monte Carlo

Um modo diferente de calcular probabilidades é por meio de experimentos computacionais. Ilustraremos esta ideia usando o chamado método de simulação de Monte Carlo. Vejamos um exemplo simples. Suponhamos que temos um caixa com 10 bolas iguais, mas com cores diferentes: 3 vermelhas, 5 azuis e 2 brancas. Vejamos com simular escolhas aleatórias, sem substituição de 5 bolas da caixa. A caixa com as bolas pode ser representada com ajuda do comando rep seguinte do R:

bolas <- rep(c("red", "blue", "white"),times = c(3,5,2))
bolas

 [1] "red"   "red"   "red"   "blue"  "blue"  "blue" 
 [7] "blue"  "blue"  "white" "white"

Selecionemos aleatoriamente, sem substituição, 5 bolas:

sample(bolas,5)

[1] "white" "blue"  "white" "blue"  "red"  

Se o comando for repetido, outra conjunto de bolas será obtido:

sample(bolas,5)

[1] "red"   "red"   "red"   "blue"  "white"

O experimento pode ser replicado quantas vezes quisermos usando o comando replicate:

B <- 6
eventos <- replicate(B, sample(bolas,5))
eventos

     [,1]   [,2]    [,3]    [,4]    [,5]   [,6]  
[1,] "blue" "red"   "blue"  "blue"  "blue" "blue"
[2,] "blue" "blue"  "red"   "blue"  "blue" "red" 
[3,] "blue" "blue"  "blue"  "white" "blue" "red" 
[4,] "red"  "blue"  "white" "blue"  "red"  "blue"
[5,] "blue" "white" "blue"  "white" "blue" "red" 

Podemos contar automaticamente o número de bolas de cada tipo obtidas da seguint forma:

tab <- table(eventos)
tab

eventos
 blue   red white 
   19     7     4 

Se selecionarmos 10 bolas o conjunto original é obtido:

sample(bolas,10)

 [1] "white" "blue"  "red"   "blue"  "white" "blue" 
 [7] "blue"  "blue"  "red"   "red"  

embora em ordem diferente. Se selecionarmos somente uma bola, é fácil entender que as probabilidades P1, P2 e P3 de obtermos bolas vermelhas, azuis e brancas devem ser dadas, P1 = 0.3, P2 = 0.5, P3 = 0.2. Façamos uma simulação de Monte Carlo com 10000 experimentos, selecionando uma bola por experimento e contando quantas são obtidas:

B <- 10000
eventos <- replicate(B,sample(bolas,1))
tab <- table(eventos)
tab

eventos
 blue   red white 
 4999  2988  2013 

Vejamos as proporções:

prop.table(tab)

eventos
  blue    red  white 
0.4999 0.2988 0.2013 

que é, aproximadamente, o resultado esperado.

Exemplo

Vejamos agora qual é a probabilidade de que tenhamos nas duas primeiras escolhas, realizadas sem substituição, obtenhamos bolas azuis. Denominemos tal probabilidade por \(Pr(BB)\). É fácil mostrar que \[ Pr(BB)= \frac{5}{10}\cdot \frac{4}{9} =\frac{2}{9} \approx 0.222 \] Vejamos como obter este resultado a partir de uma simulação de Monte Carlo. Devemos repetir o seguinte experimento diversas vezes e contar o número de vezes que obtemos “blue” “blue”:

experimento <- sample(bolas,2)
experimento

[1] "blue" "red" 

Façamos o experimento 10 vezes:

B <-10
bolasazuis <- replicate(B,{experimento <- sample(bolas,2) 
(experimento[1]=="blue" & experimento[2]=="blue")})
bolasazuis

 [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [9] FALSE  TRUE

O evento BB ocorreu 3 vezes. Calculemos a proporção de número de TRUE e número de FALSE:

mean(bolasazuis)

[1] 0.1

Façamos agora 100000 repetições:

B <-10000
bolasazuis <- replicate(B,{experimento <- sample(bolas,2) 
(experimento[1]=="blue" & experimento[2]=="blue")})
mean(bolasazuis)

[1] 0.224

de modo que reobtemos, aproximadamente, o resultado esperado.

Exemplo

Suponhamos que queremos saber agora qual é a probabilidade de obtermos os eventos \(RWB\) ao retirarmos 3 bolas, sem substituição. Sabemos que \[ Pr(RWB)=\frac{3}{10}\frac{2}{9}\frac{5}{8}=\frac{1}{24}\approx 0.04167\,. \] Vejamos como obter este resultado a partir de uma simulação de Monte Carlo, fazendo 800 000 repetições:

B <-800000
bolascoloridas <- replicate(B,{experimento <- sample(bolas,3) 
(experimento[1]=="red" & experimento[2]=="white" & experimento[3]=="blue")})
mean(bolascoloridas)

[1] 0.04179125

que é o resultado esperado.

Exemplo

Suponhamos agora que selecionamos aleatoriamente, sem substituição, 5 bolas e queremos saber qual é a probabilidade de que ao menos 2 bolas sejam vermelhas. Para resolver este problema pelo método tradicional anterior teríamos que listar muitos eventos \(RRBWW\), \(RWRRB\), etc., calcular a probabilidade de cada evento e somar. Façamos logo uma simulação de Monte Carlo, usando o comando str_count (pacote stringr) que funciona do seguinte modo:

s <- sample(bolas,5)
s

[1] "red"  "blue" "red"  "blue" "red" 

Contemos o número de bolas vermelhas na amostra:

str_count(s, "red")

[1] 1 0 1 0 1

Ou seja, quando o string é encontrado o valor 1 é atribuído à respectiva componente, caso contrário é 0. Para determinar o número de ocorrências de 1 no vetor usamos o comando:

n <- str_count(s, "red")
length(which(n==1))

[1] 3

Aplicando este comando ao nosso problema, repetindo o experimento 100 000 vezes, a proporção de casos nos quais temos duas ou mais bolas vermelhas é:

library(stringr)
B <- 100000
bolasvermelhas<-replicate(B,{experimento <- sample(bolas,5) 
n <- str_count(experimento, "red")
length(which(n==1)) >=2})
mean(bolasvermelhas)

[1] 0.49679

ou seja, aproximadamente 1/2.

Exemplo

De uma caixa com 20 bolas, sendo 6 brancas, 4 vermelhas, 7 azuis e 3 verdes são selecionadas 6, sem substituição. (i) Qual é a probabilidade de que a primeira e a última bola da seleção sejam azuis? (ii) Qual é a probabilidade de que ao menos três bolas sejam brancas?

Solução

1. Definimos a caixa com bolas:

bolas <- rep(c("red", "blue", "white", "green"),times = c(4,7,6,3))
bolas

 [1] "red"   "red"   "red"   "red"   "blue"  "blue" 
 [7] "blue"  "blue"  "blue"  "blue"  "blue"  "white"
[13] "white" "white" "white" "white" "white" "green"
[19] "green" "green"

Selecionamos 6 bolas e repetimos o experimento 900 000 vezes, contando o número de vezes em que a primeira e a sexta bola são azuis:

B <-900000
bolasazuis <- replicate(B,{experimento <- sample(bolas,6) 
(experimento[1]=="blue" & experimento[6]=="blue")})
mean(bolasazuis)

[1] 0.1102656

2. Devemos agora contabilizar todos os casos em que temos 3 ou mais bolas brancas:

B <-900000
BolasBrancas <- replicate(B,{experimento <- sample(bolas,6) 
n <- str_count(experimento, "white")
length(which(n==1)) >=3})
mean(BolasBrancas)

[1] 0.22508

Voltando ao problema do jogo blackjack, em vez de usar combinações para deduzir a exata probabilidade do 21 natural, podemos usar o método de Monte Carlo para estimar a probabilidade. Neste caso tiramos duas cartas repetidamente e contamos quantas vezes obtemos 21. Podemos usar o comando sample para retirar uma carta sem substituição:

hand <- sample(deck2, 2)
hand

[1] "Eight Diamonds" "Jack Diamonds" 

e checar se obtemos 21. Façamos este experimento 100 vezes:

B <- 100
resultados <- replicate(B,{hand <- sample(deck2,2) 
(hand[1] %in% aces & hand[2] %in% facecard) | (hand[2] %in% aces & hand[1] %in% facecard)})
resultados

  [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE
  [8] FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [15] FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [22] FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [29] FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [36] FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [43] FALSE  TRUE FALSE FALSE FALSE FALSE FALSE
 [50] FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [57] FALSE FALSE FALSE FALSE FALSE FALSE  TRUE
 [64] FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [71] FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [78] FALSE FALSE FALSE FALSE FALSE  TRUE FALSE
 [85] FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [92] FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [99] FALSE FALSE

Calculemos a proporção de número de TRUE e número de FALSE:

mean(resultados)

[1] 0.03

Façamos 100000 amostragens:

B <- 100000
resultados <- replicate(B,{hand <- sample(deck2,2) 
(hand[1] %in% aces & hand[2] %in% facecard) | (hand[2] %in% aces & hand[1] %in% facecard)})
mean(resultados)

[1] 0.04771

que é o resultado esperado.

Exemplo

Temos dois times de volei, Brasil e Argentina, que estão jogando uma séries de melhor de 7 jogos. O Brasil tem uma equipe melhor e tem 60% de chances de vencer um dado jogo.
(i) Qual é a probabilidade de que a Argentina vença ao menos um jogo?
(ii) Resolva este problema usando simulação de Monte Carlo.
Solução: Notemos que a Argentina precisa vencer ao menos um jogo dos primeiros 4 para que não seja eliminada sem vitórias. A probabililidade de que isso não aconteça é igual à probabilidade do Brasil vencer os 4 primeiros jogos, ou seja,

p_brasil_4v <- 0.6^4
p_brasil_4v

[1] 0.1296

Portanto, a probabilidade a Argentina vença ao menos um jogo entre os 4 primeiros é:

1-p_brasil_4v

[1] 0.8704

2. Façamos a simulação do experimento 100000 vezes:

B = 100000
argentina_vence <- replicate(B,{resultados <- sample(c("lose","win"), 4, replace = TRUE, prob = c(0.6,0.4)) 
any(resultados=="win")})
mean(argentina_vence)

[1] 0.8722

o que está de acordo com o cálculo exato.

O Problema dos Aniversários

Suponhamos que temos uma sala de aula com 50 alunos. Qual é a probabilidade de que dois alunos façam aniversário no mesmo dia? Usemos simulação de Monte Carlo para resolver este problema. Vamos supor que ninguém faz aniversário no dia 29 de fevereiro (o que não muda significativamente o resultado final). Aniversários podem ser representados por números de 1 a 365. Portanto, selecionemos 50 dias de 365 aleatoriamente, com substituição uma vez que dias de aniversários não são exclusivos:

n <- 50
bd <- sample(1:365,n,replace = TRUE)

Para verificar se neste conjunto de 50 dias há ao menos dois repetidos (aniversários no mesmo dia), usaremos o comando duplicated, que retorna o valor TRUE sempre que um elemento do vetor já houver aparecido antes. Por exemplo,

dup <- duplicated(c(2,3,4,2,2,5,3,2,1))
dup

[1] FALSE FALSE FALSE  TRUE  TRUE FALSE  TRUE  TRUE
[9] FALSE

Portanto, obtemos o valor TRUE para o elemento 2 na quarta e quinta componentes, onde o 2 se repete.. Além disso, obtemos TRUE para os elementos 3 e 2 na sétima e oitava componentes, onde eles repetem. Se estivermos interessados somente no caso em que a duplicação, podemos usar o comando any:

any(dup)

[1] TRUE

Apliquemos estas ideias ao problema dos aniversários, fazendo 10 amostras:

n < - 50

[1] FALSE

B <- 10
resultado <- replicate(B,{bd <- sample(1:365,n,replace = TRUE) 
any(duplicated(bd))})
resultado

 [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[10] TRUE

Notemos que em 10 experimentos obtivemos repetições em todos eles. Façamos 100 000 experimentos e calculemos a proporção de repetições:

n < - 50

[1] FALSE

B <- 100000
mesmo_dia <- replicate(B,{bd <- sample(1:365,n,replace = TRUE) 
any(duplicated(bd))})
mean(mesmo_dia)

[1] 0.97155

Portanto, em um grupo de 50 pessoas, aniversários no mesmo dia ocorrem com probabilidade de aproximadamente 97%.

Consideremos agora o problema de calcular, para um dado grupo, as chances de haver aniversário no mesmo dia com diferentes probabilidades. Criaremos uma tabela para examinar as possibilidades.

Inicialmente criamos uma função a partir da simulação de Monte Carlo definida acima:

compute_prob <- function(n,
B=100000)
{mesmo_dia <- replicate(B,{bd <- sample(1:365,n,replace = TRUE) 
any(duplicated(bd))})
mean(mesmo_dia)
}

Testemos a função

compute_prob(40)

[1] 0.89203

Definimos agora um vetor com diversos valores para \(n\):

n <- 40:60

Notemos que a nossa função não funciona como esperado:

compute_prob(n)

[1] 0.89089

Ou seja, em vez de calcular a probabilidade para todos os n a função somente usou o primeiro valor (= 40). Isso é porque nossa função espera um escalar como entrada e não um vetor. Outras funções em R normalmente aceitam vetores como argumento. Por exemplo:

sqrt(n)

 [1] 6.324555 6.403124 6.480741 6.557439 6.633250
 [6] 6.708204 6.782330 6.855655 6.928203 7.000000
[11] 7.071068 7.141428 7.211103 7.280110 7.348469
[16] 7.416198 7.483315 7.549834 7.615773 7.681146
[21] 7.745967

Devemos usar aqui o comando sapply, que permite a que a avaliação do função seja feita em todas as componentes da entrada:

sapply(n,compute_prob)

 [1] 0.89126 0.90345 0.91504 0.92419 0.93273 0.94044
 [7] 0.94757 0.95625 0.96069 0.96579 0.96950 0.97417
[13] 0.97821 0.98146 0.98411 0.98599 0.98811 0.99065
[19] 0.99162 0.99308 0.99377

Façamos um gráfico para n de 2 a 60:

n <- 2:60
prob <- sapply(n,compute_prob)
plot(n,prob)

Calculemos as probabilidades exatas em vez de usar simulações de Monte Carlo. Neste caso, é mais simples calcular a probabilidade de que o evento (aniversários no mesmo dia) não ocorra. Examinemos a probabilidade Pr(n) de que um número n de pessoas tenha aniversário numa data única. Se n=1 teremos obviamente Pr(1)=1. A probabilidade de que a segunda pessoa tenha aniversário único, dado que a primeira pessoa tinha aniversário único é 363/365. A probabilidade de que a terceira pessoa tenha aniversário único, dado que as duas primeiras tinham aniversário único é 363/365. Assim, as probabilidades 2 e 3 pessoas tenham aniversário único são dada por \[ Pr(2) = 1 \times \frac{364}{365}\,, \quad Pr(3) = 1 \times \frac{364}{365} \times \frac{363}{365}\,. \] De modo geral, \[ Pr(n) = 1 \times \frac{364}{365} \times \frac{363}{365} \times \cdots \times \frac{365-n+1}{365}\,. \] Façamos a implementação em R:

p_exata <- function(n){
  p_unico <- seq(365, 365-n+1)/365
  1 - prod(p_unico)
}

Façamos um teste:

p_exata(40)

[1] 0.8912318

Podemos agora avaliar esta função nos vários valores de n:

eprob <- sapply(n,p_exata)
plot(n, eprob)

Comparemos este resultado gráfico com aquele obtido com o método de Monte Carlo, usando um curva contínua para o resultado exato:

plot(n, prob)
lines(n, eprob, col="blue")

Tal resultado mostra uma boa concordância entre os métodos.

Número de Amostras

Uma questão crucial em simulações de Monte Carlo envolve a escolha do número de experimentos a serem feitos. Definimos um intervalo de variação de B:

B <- 10^seq(1, 5, len = 100)

e agora uma função que calcula as probabilidades para cada simulação. Fixemos n=30.

compute_prob <- function(B, n = 30){
  same_day <- replicate (B, {
    bdays <- sample(1:365, n, replace=TRUE)
    any(duplicated(bdays))
  })
  mean(same_day)
}

Façamos agora a avaliação da função para todos os valores de B e plotemos o resultado:

prob <- sapply(B, compute_prob)
plot(log10(B), prob, type ="l")

Notemos que a partir de \(B = 10^3\) os resultados começam a estabilizar.

Exercícios

1. Retiramos, com substituição, 5 bolas de uma caixa que contém 4 bolas brancas, 6 azuis e 7 amarelas. Supondo que as 5 primeiras bolas foram todas amarelas, qual é a probabilidade de que ao retirar a sexta bola encontremos novamente uma bola amarela? R: 0.47

2. Jogamos um dado honesto 6 vezes. Qual é a probabilidade de não sair um 6 em nenhuma jogada? R: 0.33

3. Um time de futebol precisa fazer 7 pontos em 5 jogos para não ser rebaixado para a liga inferior. Nos jogos que fará estima-se que ele tem 10% de chances de vencer, 50% de empatar e 40% de perder. Qual é a probabilidade do time não ser rebaixado. Calcule exatamente e usando uma simulação de Monte Carlo.