Probabilidade Discreta e Simulações de Monte Carlo

Probabilidade Discreta

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íveis 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.

Definindo cartas de um baralho em R

Nesta seção 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"   

Exemplo. Retirando uma carta do baralho

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. Confirmaremos este resultado em R.

Permutações em R

O comando permutations() 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, como veremos a seguir.

Exemplo. Permutação básica

façamos uma permutação (ou arranjo) de 4 números, 3,4,5,6 em grupos de 3 elementos. De modo mais concreto, se temos uma caixa com 4 bolas marcadas com estes números, quais são os possíveis resultados ao selecionarmos 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] 10 13 22 14  4 21

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

m[indice,]

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

Exemplo. Permutação

Simule a seleção aleatória de de 5 números (0 a 9) de 7 dígitos.
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] 187215 360073 178116 292975 259809

numeros[indice,]

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

Exemplo. Probabilidade condicional

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)} \]

Combinações em R

Suponhamos agora que a ordem não seja importante. 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

Exemplo. Blackjack

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. 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ção 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. 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] "blue"  "blue"  "red"   "white" "blue" 

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

sample(bolas,5)

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

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,] "white" "red"  "blue"  "red"   "blue"  "blue"
[2,] "blue"  "blue" "red"   "white" "red"   "blue"
[3,] "blue"  "blue" "red"   "red"   "blue"  "red" 
[4,] "white" "red"  "red"   "blue"  "white" "red" 
[5,] "red"   "blue" "white" "blue"  "white" "red" 

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

tab <- table(eventos)
tab

eventos
 blue   red white 
   12    12     6 

Se selecionarmos 10 bolas o conjunto original é obtido:

sample(bolas,10)

 [1] "blue"  "white" "red"   "blue"  "red"   "blue" 
 [7] "blue"  "white" "blue"  "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 
 5110  2979  1911 

Vejamos as proporções:

prop.table(tab)

eventos
  blue    red  white 
0.5110 0.2979 0.1911 

que é, aproximadamente, o resultado esperado.

Exemplo. Simulação de Monte Carlo e condição lógica (1)

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] "white" "white"

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  TRUE
 [9] FALSE FALSE

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.226

de modo que reobtemos, aproximadamente, o resultado esperado.

Exemplo. Simulação de Monte Carlo e condição lógica (2)

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.041495

que é o resultado esperado.

Exemplo. Simulação de Monte Carlo e condição lógica (3)

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] "blue"  "blue"  "blue"  "blue"  "white"

Contemos o número de bolas vermelhas na amostra:

str_count(s, "red")

[1] 0 0 0 0 0

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] 0

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.50112

ou seja, aproximadamente 1/2.

Exemplo. Simulação de Monte Carlo e condição lógica (4)

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.1106511

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.2254611

Exemplo. Blackjack

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] "King Clubs" "Jack Clubs"

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]  TRUE FALSE FALSE FALSE FALSE FALSE FALSE
 [50]  TRUE FALSE FALSE FALSE FALSE FALSE FALSE
 [57] FALSE FALSE FALSE FALSE FALSE  TRUE FALSE
 [64] FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [71] FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [78] FALSE FALSE FALSE  TRUE FALSE FALSE 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.04

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.04845

que é o resultado esperado.

Exemplo. Jogo de volei

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.86982

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

Exemplo. Outra série de jogos de volei

Dois times de volei, Brasil e Sérvia, estão jogando uma série de 7 jogos. O primeiro a vencer 4 jogos vence a série. Ambos os times são equilibrados e ambos têm iguais chances de vencer cada jogo. Se o Brasil vencer o primeiro jogo, qual é a chance de que ele vença a série? Calcule exatamente e depois faça uma simulação de Monte Carlo.

Solução

Calculemos exatamente. O número de jogos restantes é n = 6. Definimos uma variável chamada resultados, que é o vetor de possíveis resultados, sendo 0 uma derrota e 1 uma vitória

n <-6
resultados <- c(0,1)

Criamos uma lista l com todos os possíveis resultados dos jogos remanescentes:

l <- replicate(n,list(resultados))

Criamos um data frame named possibilidades que contém todas as combinações de todos os possíveis resultados para os jogos remanescentes:

possibilidades <- expand.grid(l)
possibilidades

Criamos um vetor chamado resultados que indica se cada linha no data frame possibilidades contém o número suficiente de vitórias para que o time do Brasil vença a série:

resultados <- rowSums(possibilidades)>3
resultados

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

Calculamos em resultados a proporção de casos (TRUE) em que o Brasil vence a série:

mean(resultados)

[1] 0.34375

Façamos agora uma simulação de Monte Carlo. Criamos um objeto chamado resultados que replica para B = 10000 iterações uma série simulada e determina se a série contém ao menos 4 vitórias, o que daria vitória ao time do Brasil e depois calculemos a proporção de casos TRUE

B <- 10000
resultados <- replicate(B,sum(replicate(6, sample(c(0,1),1)))>3)
mean(resultados)

[1] 0.3435

Deste modo, confirmamos o resultado obtido exatamente.

Exemplo. Séries de jogos com diferentes probabilidades

Dois times, A e B jogam uma série de 7 jogos. A probabilidade de A vencer B é p. Determine as probabilidades de que A vença a série para p = 0.5, 0.525, 0.550, … 0.95, e apresente o resultado em um gráfico.

Solução

Definimos o vetor de probabilidades:

p <- seq(0.5, 0.95, 0.025)

Definimos a função:

prob_win <- function(p){
  B <- 10000
  result <- replicate(B, {
    b_win <- sample(c(1,0), 7, replace = TRUE, prob = c(1-p, p))
    sum(b_win)>=4
    })
    mean(result)
}

Façamos agora a avaliação da função para cada valor de p e o correspondente gráfico:

Pr <- sapply(p,prob_win)
plot(p,Pr)

Exemplo. Diferentes números de jogos em cada série

Consideremos um problema semelhante ao anterior, mas agora com probabilidade de que A vença qualquer jogo fixa em p=0.75, mas com diferentes séries com diferentes números de jogos.

Solução.

Devemos agora deixar o número de jogos N como argumento:

prob_win <- function(N, p=0.75){
      B <- 10000
      result <- replicate(B, {
        b_win <- sample(c(1,0), N, replace = TRUE, prob = c(1-p, p))
        sum(b_win)>=(N+1)/2
        })
            mean(result)
    }

Definimos os possíveis valores para N como números ímpares de 1 a 25:

N <- seq(1, 25, 2)
N

 [1]  1  3  5  7  9 11 13 15 17 19 21 23 25

Façamos agora a avaliação da função para cada valor de N e o correspondente gráfico:

Pr <- sapply(N,prob_win)
plot(N,Pr)

Exemplo. 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.96977

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

Exemplo. O problema dos aniversários com diversos grupos

Consideremos agora o problema de calcular, para um dado grupo, as chances de haver aniversário no mesmo dia para diferentes grupos com diferente número de pessoas. 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.8909

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.89049

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.89106 0.90290 0.91373 0.92242 0.93318 0.94118
 [7] 0.94794 0.95505 0.96021 0.96646 0.97077 0.97439
[13] 0.97704 0.98124 0.98390 0.98626 0.98817 0.99028
[19] 0.99119 0.99267 0.99445

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 nas simulações de Monte Carlo

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.

Exemplo. O Problema de Monty Hall

O problema de Monty Hall pode ser formulado da seguinte maneira:

Suponhamos que você está num concurso e tem que escolher entre 3 portas: atrás de uma das portas está um carro, e atrás das outras bodes. Você escolhe uma porta, digamos a porta 1 e o apresentador, que sabe onde está o carro, vai abrir outra porta, 3, por exemplo, onde há um bode. Ele então lhe pergunta: “você quer ficar com a porta 1 ou quer trocar pela porta 2?” Há alguma vantagem em se trocar de porta? Abaixo está o código referente à escolha da primeira porta.

B <- 10000 #número de experimentos
mesma_porta <- replicate(B, {
portas <- as.character(1:3) #define o vetor com components "1", "2", "3". 
premio <- sample(c("carro", "bode", "bode")) #mistura aleatoriament ordem do string "carro", "bode", "bode"
porta_premiada <- portas[premio == "carro"] #determina o índice da porta premiada
minha_escolha<- sample(portas, 1) #faz uma escolha aleatória entre "1", "2", "3"
#show <- sample(portas[!portas %in% c(minha_escolha, porta_premiada)],1)
mesma_porta <- minha_escolha #mantém a escolha original
mesma_porta == porta_premiada #verifica se a escolha da porta corresponde ao prêmio (TRUE ou FALSE)
})
mean(mesma_porta) #calcule a proporção de valores TRUE

[1] 0.3328

Este resultado indica que a probabilidade é 1/3, o que é esperado, pois nada mudou relativamente à escolha original de 1 porta entre 3. Vejamos agora o que ocorre se escolhemos a estratégia de trocar a escolha:

B <- 10000 #número de experimentos
troca_porta <- replicate(B, {
portas <- as.character(1:3) #define o vetor com components "1", "2", "3". 
premio <- sample(c("carro", "bode", "bode")) #mistura aleatoriament ordem do string "carro", "bode", "bode"
porta_premiada <- portas[premio == "carro"] #determina o índice da porta premiada
minha_escolha<- sample(portas, 1) #faz uma escolha aleatória entre "1", "2", "3"
show <- sample(portas[!portas %in% c(minha_escolha, porta_premiada)],1) #mostra a porta diferente da originalmente escolhida e que não contém o carro
troca_porta <- portas[!portas%in%c(minha_escolha, show)]#escolhe a porta que não é a mostrada (pois não contém prêmio e também não é a original)
troca_porta == porta_premiada #verifica se a escolha da porta corresponde ao prêmio (TRUE ou FALSE)
})
mean(troca_porta) #calcule a proporção de valores TRUE

[1] 0.6686

Ou seja, a probabilidade agora é 2/3. Isso é esperado, pois temos agora a informação de qual é uma das portas que não contém o prêmio.

Exercícios

Em todos os problemas abaixo, calcule exatamente e simule com Monte Carlo.

1. Uma caixa contém 15 bolas brancas e 25 bolas azuis. Duas bolas são selecionadas aleatóriamente, sem substituição. Seja A o evento correspondente ao fato da primeira bola ser branca e B o evento associado ao fato da segunda bola ser azul. Determine \(Pr(A)\), \(Pr(B|A)\), \(Pr(A \cup B)\), \(Pr(A \cap B)\).

2. (MR p.61 [1]) O circuito mostrado na figura abaixo funciona se e somente se há um caminho de dispositivos funcionais da esquerda para a direita. Suponha que todos os dispositivos falham independentemente que a probabilidade de falha em cada dispositivo é indicada pelo número em cada caixa. Qual é a probabilidade de que o circuito opere?

3. 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

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

5. Jogamos 2 dados honestos. Qual é a probabilidade de obtermos mais de 10 pontos em 4 jogadas?

6. Resolva os problemas anteriores supondo que os dados são viciados do seguinte modo: o lado 6 tem probabilidade 1.1/6 e os outros lados têm probabilidades iguais.

7. 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.

8. Em um baralho completo, ao retirarmos 5 cartas, determine a probabilidade conseguirmos cada um dos seguintes resultados:

1. Dois pares
   
2. Trinca
   
3. Sequência
   
4. Flush
   
5. Full House
   
6. Quadra
   
7. Straigh Flash
   
8. Royal Flash.

Referências

1. Douglas C. Montgomery e George C. Runger. Applied Statistics and Probability for Engineers, 6th ed., Wiley, 2014