Introdução à Simulação

Fernando Deeke Sasse
CCT - UDESC

1. Permutação e computação de probabilidades por meio de contagem

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

Error in permutations(4, 3, c(3, 4, 5, 6)) : 
  could not find function "permutations"

Notemos que uma ou mais linhas específicas de m podem ser obtidas da seguinte maneira. Por exemplo, selecionemos as linhas 12,13 e 14:

Podemos também usar um vetor de strings:

m <- permutations(4,3,c("branca","preta","azul","vermelha"))
m

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

Como vemos da matriz de saída em qualquer um dos casos 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 5 amostragens. Para isso inicialmente escolhemos aleatóriamente 4 números entre os n=24:

indice <- sample(n,4)
indice

[1]  5  4 19  6

Cada uma destas cinco escolhas corresponde a um resultado de um experimento:

m[indice,]

     [,1]       [,2]       [,3]      
[1,] "azul"     "vermelha" "branca"  
[2,] "azul"     "preta"    "vermelha"
[3,] "vermelha" "azul"     "branca"  
[4,] "azul"     "vermelha" "preta"   

A probabilidade de numa amostragem obtermos a bola azul é obviamente 1/4 que corresponde a 6/24 no experimento que consiste em coletar 3 bolas aleatoriamente. S

Suponhamos adicionamos agora mais uma bola branca ao conjunto original de 4 bolas. Verifiquemos agora as possibilidades correspondentes à seleção de 3 bolas:

bolas = c("branca1","preta","azul","vermelha", "branca2")
m <- permutations(5,3,bolas)
m

      [,1]       [,2]       [,3]      
 [1,] "azul"     "branca1"  "branca2" 
 [2,] "azul"     "branca1"  "preta"   
 [3,] "azul"     "branca1"  "vermelha"
 [4,] "azul"     "branca2"  "branca1" 
 [5,] "azul"     "branca2"  "preta"   
 [6,] "azul"     "branca2"  "vermelha"
 [7,] "azul"     "preta"    "branca1" 
 [8,] "azul"     "preta"    "branca2" 
 [9,] "azul"     "preta"    "vermelha"
[10,] "azul"     "vermelha" "branca1" 
[11,] "azul"     "vermelha" "branca2" 
[12,] "azul"     "vermelha" "preta"   
[13,] "branca1"  "azul"     "branca2" 
[14,] "branca1"  "azul"     "preta"   
[15,] "branca1"  "azul"     "vermelha"
[16,] "branca1"  "branca2"  "azul"    
[17,] "branca1"  "branca2"  "preta"   
[18,] "branca1"  "branca2"  "vermelha"
[19,] "branca1"  "preta"    "azul"    
[20,] "branca1"  "preta"    "branca2" 
[21,] "branca1"  "preta"    "vermelha"
[22,] "branca1"  "vermelha" "azul"    
[23,] "branca1"  "vermelha" "branca2" 
[24,] "branca1"  "vermelha" "preta"   
[25,] "branca2"  "azul"     "branca1" 
[26,] "branca2"  "azul"     "preta"   
[27,] "branca2"  "azul"     "vermelha"
[28,] "branca2"  "branca1"  "azul"    
[29,] "branca2"  "branca1"  "preta"   
[30,] "branca2"  "branca1"  "vermelha"
[31,] "branca2"  "preta"    "azul"    
[32,] "branca2"  "preta"    "branca1" 
[33,] "branca2"  "preta"    "vermelha"
[34,] "branca2"  "vermelha" "azul"    
[35,] "branca2"  "vermelha" "branca1" 
[36,] "branca2"  "vermelha" "preta"   
[37,] "preta"    "azul"     "branca1" 
[38,] "preta"    "azul"     "branca2" 
[39,] "preta"    "azul"     "vermelha"
[40,] "preta"    "branca1"  "azul"    
[41,] "preta"    "branca1"  "branca2" 
[42,] "preta"    "branca1"  "vermelha"
[43,] "preta"    "branca2"  "azul"    
[44,] "preta"    "branca2"  "branca1" 
[45,] "preta"    "branca2"  "vermelha"
[46,] "preta"    "vermelha" "azul"    
[47,] "preta"    "vermelha" "branca1" 
[48,] "preta"    "vermelha" "branca2" 
[49,] "vermelha" "azul"     "branca1" 
[50,] "vermelha" "azul"     "branca2" 
[51,] "vermelha" "azul"     "preta"   
[52,] "vermelha" "branca1"  "azul"    
[53,] "vermelha" "branca1"  "branca2" 
[54,] "vermelha" "branca1"  "preta"   
[55,] "vermelha" "branca2"  "azul"    
[56,] "vermelha" "branca2"  "branca1" 
[57,] "vermelha" "branca2"  "preta"   
[58,] "vermelha" "preta"    "azul"    
[59,] "vermelha" "preta"    "branca1" 
[60,] "vermelha" "preta"    "branca2" 

Temos agora 60 possibilidades:

n <- nrow(m)
n

[1] 60

Definimos as colunas:

primeira_bola <- m[,1]
segunda_bola <- m[,2]
terceira_bola <- m[,3]

1. Em quais casos a primeira bola é azul?

h<-primeira_bola %in% "azul"
h

 [1]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE FALSE FALSE
[16] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[31] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[46] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE

Obviamente isso acontece 12 vezes:

sum(h)

[1] 12

2. Em quais casos a primeira bola é branca? Note que a primeira bola é branca quando selecionamos a bola branca1 ou a bola branca2, de modo que o número deve ser 24. Para sistematizar isso definimos o vetor:

brancas<- c("branca1", "branca2")

Portanto,

h<-primeira_bola %in% brancas
h

 [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE  TRUE  TRUE
[16]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
[31]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[46] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE

Obviamente isso acontece 24 vezes:

sum(h)

[1] 24

3. Qual é a probabilidade de que a primeira bola seja branca ou azul?

N<-sum((primeira_bola %in% brancas) | primeira_bola %in% "azul")
N

[1] 36

A probabilidade é dada por

p=N/n
p

[1] 0.6

4. Qual é a probabilidade de que a primeira bola seja branca e a terceira azul?

N<-sum((primeira_bola %in% brancas) & terceira_bola %in% "azul")
N

[1] 6

A probabilidade é dada por

p=N/n
p

[1] 0.1

Exercício 1.1 Em uma caixa temos 2 bolas brancas, 3 azuis, uma preta, uma vermelha e uma amarela. Qual é a probabilidade de que em uma amostragem na qual selecionamos 4 bolas obtenhamos a primeira bola amarela, a segunda branca, a terceira e quarta azuis?

Solução. Definimos nosso conjunto de bolas e exibimos o conjunto de todas as possíveis permutações (arranjos):

bolas = c("branca1", "branca2","preta","azul1","azul2", "azul3", "vermelha", "amarela")
m <- permutations(8,4,bolas)
n <- nrow(m)
n

[1] 1680

Temos portanto 1680 possibilidades. Mostramos abaixo seis delas:

head(m)

     [,1]      [,2]    [,3]    [,4]      
[1,] "amarela" "azul1" "azul2" "azul3"   
[2,] "amarela" "azul1" "azul2" "branca1" 
[3,] "amarela" "azul1" "azul2" "branca2" 
[4,] "amarela" "azul1" "azul2" "preta"   
[5,] "amarela" "azul1" "azul2" "vermelha"
[6,] "amarela" "azul1" "azul3" "azul2"   

Definimos as colunas:

primeira_bola <- m[,1]
segunda_bola <- m[,2]
terceira_bola <- m[,3]
quarta_bola <- m[,4]

Vejamos quantos eventos resultam em “primeira bola amarela, a segunda branca, a terceira ou quarta azuis”. Definimos antes os conjuntos de bolas brancas e azuis:

brancas<- c("branca1", "branca2")
azuis<- c("azul1","azul2", "azul3")

Portanto,

N<-sum((primeira_bola %in% "amarela") & (segunda_bola %in% brancas) & (terceira_bola %in% azuis | quarta_bola %in% azuis))      
N

[1] 48

Portanto, a probabilidade do evento especificado ocorrer é dada por

p <- N/n
p

[1] 0.02857143

Exercício 1.2 Em uma caixa temos 2 bolas brancas, 3 azuis, uma preta, uma vermelha e uma amarela. Qual é a probabilidade de que em uma amostragem na qual selecionamos 2 bolas obtenhamos as duas primeiras bolas azuis?

Solução. Definimos nosso conjunto de bolas e exibimos o conjunto de todas as possíveis permutações (arranjos):

bolas = c("branca1", "branca2","preta","azul1","azul2", "azul3", "vermelha", "amarela")
m <- permutations(8,2,bolas)
n <- nrow(m)
n

[1] 56

Temos portanto 56 possibilidades. Definimos as colunas:

primeira_bola <- m[,1]
segunda_bola <- m[,2]

Vejamos quantos eventos resultam em “primeira bola azul e segunda azul”:

azuis<- c("azul1","azul2", "azul3")

Portanto,

N<-sum((primeira_bola %in% azuis) & (segunda_bola %in% azuis))      
N

[1] 6

Portanto, a probabilidade do evento especificado ocorrer é dada por

p <- N/n
p

[1] 0.1071429

2. Combinações

Quando num processo de amostragem, no caso dos exemplos acima, a ordem dos eventos não é importante, podemos usar combinações. 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 2.1 Suponhamos que uma caixa temos 2 bolas brancas, 3 azuis, uma preta, uma vermelha e uma amarela. Qual é a probabilidade de que numa seleção de 4 bolas obtenhamos 2 bolas azuis e uma bola amarela?

Solução. Definimos nosso conjunto de bolas e exibimos o conjunto de todas as possíveis combinações:

bolas = c("branca1", "branca2","preta","azul1","azul2", "azul3", "vermelha", "amarela")
m <- combinations(8,4,bolas)
m

      [,1]      [,2]      [,3]      [,4]      
 [1,] "amarela" "azul1"   "azul2"   "azul3"   
 [2,] "amarela" "azul1"   "azul2"   "branca1" 
 [3,] "amarela" "azul1"   "azul2"   "branca2" 
 [4,] "amarela" "azul1"   "azul2"   "preta"   
 [5,] "amarela" "azul1"   "azul2"   "vermelha"
 [6,] "amarela" "azul1"   "azul3"   "branca1" 
 [7,] "amarela" "azul1"   "azul3"   "branca2" 
 [8,] "amarela" "azul1"   "azul3"   "preta"   
 [9,] "amarela" "azul1"   "azul3"   "vermelha"
[10,] "amarela" "azul1"   "branca1" "branca2" 
[11,] "amarela" "azul1"   "branca1" "preta"   
[12,] "amarela" "azul1"   "branca1" "vermelha"
[13,] "amarela" "azul1"   "branca2" "preta"   
[14,] "amarela" "azul1"   "branca2" "vermelha"
[15,] "amarela" "azul1"   "preta"   "vermelha"
[16,] "amarela" "azul2"   "azul3"   "branca1" 
[17,] "amarela" "azul2"   "azul3"   "branca2" 
[18,] "amarela" "azul2"   "azul3"   "preta"   
[19,] "amarela" "azul2"   "azul3"   "vermelha"
[20,] "amarela" "azul2"   "branca1" "branca2" 
[21,] "amarela" "azul2"   "branca1" "preta"   
[22,] "amarela" "azul2"   "branca1" "vermelha"
[23,] "amarela" "azul2"   "branca2" "preta"   
[24,] "amarela" "azul2"   "branca2" "vermelha"
[25,] "amarela" "azul2"   "preta"   "vermelha"
[26,] "amarela" "azul3"   "branca1" "branca2" 
[27,] "amarela" "azul3"   "branca1" "preta"   
[28,] "amarela" "azul3"   "branca1" "vermelha"
[29,] "amarela" "azul3"   "branca2" "preta"   
[30,] "amarela" "azul3"   "branca2" "vermelha"
[31,] "amarela" "azul3"   "preta"   "vermelha"
[32,] "amarela" "branca1" "branca2" "preta"   
[33,] "amarela" "branca1" "branca2" "vermelha"
[34,] "amarela" "branca1" "preta"   "vermelha"
[35,] "amarela" "branca2" "preta"   "vermelha"
[36,] "azul1"   "azul2"   "azul3"   "branca1" 
[37,] "azul1"   "azul2"   "azul3"   "branca2" 
[38,] "azul1"   "azul2"   "azul3"   "preta"   
[39,] "azul1"   "azul2"   "azul3"   "vermelha"
[40,] "azul1"   "azul2"   "branca1" "branca2" 
[41,] "azul1"   "azul2"   "branca1" "preta"   
[42,] "azul1"   "azul2"   "branca1" "vermelha"
[43,] "azul1"   "azul2"   "branca2" "preta"   
[44,] "azul1"   "azul2"   "branca2" "vermelha"
[45,] "azul1"   "azul2"   "preta"   "vermelha"
[46,] "azul1"   "azul3"   "branca1" "branca2" 
[47,] "azul1"   "azul3"   "branca1" "preta"   
[48,] "azul1"   "azul3"   "branca1" "vermelha"
[49,] "azul1"   "azul3"   "branca2" "preta"   
[50,] "azul1"   "azul3"   "branca2" "vermelha"
[51,] "azul1"   "azul3"   "preta"   "vermelha"
[52,] "azul1"   "branca1" "branca2" "preta"   
[53,] "azul1"   "branca1" "branca2" "vermelha"
[54,] "azul1"   "branca1" "preta"   "vermelha"
[55,] "azul1"   "branca2" "preta"   "vermelha"
[56,] "azul2"   "azul3"   "branca1" "branca2" 
[57,] "azul2"   "azul3"   "branca1" "preta"   
[58,] "azul2"   "azul3"   "branca1" "vermelha"
[59,] "azul2"   "azul3"   "branca2" "preta"   
[60,] "azul2"   "azul3"   "branca2" "vermelha"
[61,] "azul2"   "azul3"   "preta"   "vermelha"
[62,] "azul2"   "branca1" "branca2" "preta"   
[63,] "azul2"   "branca1" "branca2" "vermelha"
[64,] "azul2"   "branca1" "preta"   "vermelha"
[65,] "azul2"   "branca2" "preta"   "vermelha"
[66,] "azul3"   "branca1" "branca2" "preta"   
[67,] "azul3"   "branca1" "branca2" "vermelha"
[68,] "azul3"   "branca1" "preta"   "vermelha"
[69,] "azul3"   "branca2" "preta"   "vermelha"
[70,] "branca1" "branca2" "preta"   "vermelha"

Todas as 70 possibilidades estão listadas acima. De fato,

n <- nrow(m)
n

[1] 70

Definimos as colunas:

primeira_bola <- m[,1]
segunda_bola <- m[,2]
terceira_bola <- m[,3]
quarta_bola <- m[,4]

No caso de permutações tínhamos um espaço amostral de 1680 eventos. Vejamos quantos eventos resultam em “duas azuis”.

azuis <- c("azul1", "azul2", "azul3")

Para cada linha de um a n verificamos se o número de ocorrências de azul é igual a 2:

for(i in 1:n){
   a[i] <- sum(m[i,] %in% azuis) == 2 
}
print(a[1:n])

 [1] FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE
[16]  TRUE  TRUE  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[31] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
[46]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE FALSE FALSE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE
[61]  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE

Portanto o número de vezes em que temos 2 bolas azuis é dado por

N<-sum(a[1:n])
N

[1] 30

A probabilidade deste evento é dada por

p=N/n
p

[1] 0.4285714

Exemplo 2.2

Qual é a probabilidade no problema anterior, de encontrarmos duas bolas azuis e uma amarela?

Para cada linha de um a n verificamos se o número de ocorrências de azul é igual a 2 e o de amarelas é igual a 1:

for(i in 1:n){
   a[i] <- sum(m[i,] %in% azuis) == 2 & sum(m[i,] %in% "amarela")==1
}
print(a[1:n])

 [1] FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE
[16]  TRUE  TRUE  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[31] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[46] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[61] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE

Portanto o número de vezes em que temos 2 bolas azuis e uma amarela é dado por

N<-sum(a[1:n])
N

[1] 12

A probabilidade deste evento é dada por

p=N/n
p

[1] 0.1714286

3. Simulação de Monte Carlo

3.1 Introdução

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. Para uma abordagem mais detalhada deste assunto veja, por exemplo as referências [1,2,3,4].

Um modelo utilizado frequentemente para ilustrar os conceitos de inferência estatística é o do caixa. Um exemplo simples de variável aleatória é aquele que resulta quando retiramos um ticket de uma caixa (sem olhar ou aleatoriamente). Por exemplo, o lançamento de uma moeda e a determinação do lado que resulta para cima, cara (H) ou coroa (T), pode ser representado pela retirada de um ticket de uma caixa que contém dois tickets, um marcado com H e outro com T. Denotemos a caixa pelo seu conteúdo, [H,T].

Se aleatoriamente apanhamos um ticket da caixa [H,T] o ticket retirado pode ser H ou T com igual probabilidade. Suponhamos que o ticket retirado resultou ser H. Se não devolvemos o ticket original então a caixa será agora [ ,T] e se tirarmos outro ticket da caixa, então ele será T, com certeza. Se devolvermos o ticket H à caixa, então a caixa voltará a ser [H,T] novamente e uma nova retirda poderá resultar em H ou T, com igual probabilidade. Em outras palavras, se não fizermos substituição, a segunda retirada será dependente da primeira. Se fizermos substituição, a segunda retirada será independente da primeira.

O processo de substituição correponde a termos uma caixa com um número muito grande de tickets marcados H e T (em número igual). Se fizermos duas retiradas, a segunda será praticamente independente da primeira.

Criaremos tal caixa em R a seguir:

Caixa <- c("H", "T")
Caixa

[1] "H" "T"

Simulemos a retirada de 5 tickets vezes da caixa, com subtituição. Equivalentemente, poderíamos imaginar o lançamento de uma moeda 5 vezes para observar o resultado cara (H) ou coroa (T).

# Retirada de tickets da caixa
sample(x = Caixa, size = 5, replace = TRUE)

[1] "H" "H" "T" "T" "H"

Cada vez que o comando acima for executado um resultado diferente será obtido.

Os argumentos da função sample são: 1. O objeto a ser amostrado; 2. O tamanho da amostra; 3. Se a amostragem é com substituição ou não.

Ao usarmos replace = TRUE estamos supondo que todas as amostragens são independentes. Isso é equivalente à amostragem de uma população muito grande. Se replace = False (default) as amostragens não são independentes. Neste caso, se o tamanho da amostragem for maior do o número de elementos na caixa teremos um erro.

Embora cada vez que o comando sample seja executado uma diferente sequência de resultados seja gerada, em algumas situações, desejamos que uma dada simulação seja reprodutível, ou seja, que as diferentes sequêncais de resultados sempre sejam as mesmas. Para isso procedemos do seguinto modo:

set.seed(30) # para reproducibilidade
sample(Caixa, 5, replace = TRUE)

[1] "H" "H" "H" "H" "H"

Cada vez que o bloco acima for executado o resultado será o mesmo. O bloco abaixo dará um resultado diferente, mas que será sempre este, caso seja executado do bloco acima.

sample(Caixa, 5, replace = TRUE)

[1] "T" "H" "T" "T" "H"

Em uma simulação estabelecemos as regras de uma processo aleatório. Os números aleatórios gerados no computador resultam numa saída de acordo com tais regras.

Exemplo 3.1

Simulemos o experimento que consiste em retirar de uma caixa uma bola que contém bolas pretas e brancas.

O vetor caixa representa a caixa com as duas bolas, uma preta, outra branca. A função sample() realiza o experimento que consiste em retirar aleatóriamente uma bola da caixa e informar se o resultado é preta ou branca. A opção size=1 indica que só um experimento é realizado. A opção replace = TRUE aqui não faz diferença no caso de um experimento, mas é fundamental se fizermos mais experimentos com o mesmo comando, pois ele nos diz que após um dado experimento a bola selecionada é colocada de volta na caixa. Cada vez que o a função sample é chamada um dos resultado pode aparecer com probabilidade 1/2. Em vez de executarmos o comando diversas vezes podemos simplesmente especificar a quantidade de experimentos em size. Realizemos o experimento 100 vezes:

caixa <- c("preta", "branca")
resultados <- sample(caixa, size = 100, replace = TRUE)
resultados

  [1] "preta"  "branca" "preta"  "preta"  "preta"  "preta"  "preta"  "preta"  "branca" "preta" 
 [11] "branca" "preta"  "preta"  "branca" "preta"  "branca" "preta"  "preta"  "preta"  "preta" 
 [21] "preta"  "branca" "branca" "preta"  "preta"  "branca" "branca" "preta"  "branca" "preta" 
 [31] "preta"  "preta"  "preta"  "preta"  "preta"  "branca" "preta"  "preta"  "branca" "branca"
 [41] "preta"  "preta"  "branca" "branca" "branca" "branca" "preta"  "preta"  "branca" "branca"
 [51] "branca" "preta"  "branca" "branca" "branca" "preta"  "branca" "branca" "preta"  "preta" 
 [61] "preta"  "preta"  "branca" "branca" "preta"  "preta"  "preta"  "preta"  "branca" "branca"
 [71] "preta"  "preta"  "branca" "branca" "preta"  "branca" "preta"  "preta"  "preta"  "preta" 
 [81] "preta"  "branca" "branca" "branca" "preta"  "preta"  "preta"  "branca" "branca" "branca"
 [91] "preta"  "preta"  "branca" "branca" "preta"  "branca" "preta"  "preta"  "branca" "branca"

Para computar o número de bolas pretas e brancas resultantes usamos a função table():

table(resultados)

resultados
branca  preta 
    43     57 

Nos resultados acima fizemos a suposição tácita de que a probabilidade de que seja selecionada coda uma das bolas é a mesma. Este problema poderia ter sido formulado também em termos do lançamento de uma moeda honesta (que ao ser lançada tem probabilidade de resultar cara ou coroa com propabilidade 1/2). Reformulando este problema agora como o do lançamento de uma moeda, suponhamos que moeda seja “desonesta”, ou seja, que tenha probabilidade 60% de resultar cara e 40% de resultar em coroa. Basta que adicionemos o argumento prob à função sample:

moeda <- c("cara", "coroa")
resultados <- sample(moeda, size = 100, replace = TRUE, prob=c(0.6,0.4))
resultados

  [1] "coroa" "cara"  "coroa" "cara"  "coroa" "cara"  "cara"  "cara"  "coroa" "cara"  "coroa"
 [12] "cara"  "cara"  "cara"  "coroa" "cara"  "coroa" "coroa" "coroa" "cara"  "coroa" "coroa"
 [23] "cara"  "cara"  "coroa" "cara"  "cara"  "cara"  "cara"  "cara"  "coroa" "cara"  "cara" 
 [34] "cara"  "coroa" "cara"  "cara"  "coroa" "cara"  "cara"  "coroa" "coroa" "cara"  "cara" 
 [45] "cara"  "coroa" "cara"  "coroa" "coroa" "coroa" "cara"  "coroa" "cara"  "coroa" "coroa"
 [56] "coroa" "cara"  "cara"  "coroa" "cara"  "coroa" "coroa" "cara"  "coroa" "cara"  "coroa"
 [67] "coroa" "cara"  "cara"  "cara"  "coroa" "coroa" "cara"  "cara"  "coroa" "cara"  "cara" 
 [78] "cara"  "cara"  "coroa" "coroa" "cara"  "cara"  "coroa" "coroa" "coroa" "cara"  "cara" 
 [89] "cara"  "coroa" "cara"  "coroa" "cara"  "cara"  "coroa" "coroa" "cara"  "cara"  "cara" 
[100] "cara" 

Contemos os resultados:

table(resultados)

resultados
 cara coroa 
   57    43 

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

B <- 5
eventos <- replicate(B, sample(moeda, size = 100, replace = TRUE, prob=c(0.6,0.4)))
summary(eventos)

     V1         V2         V3         V4         V5    
 cara :63   cara :68   cara :63   cara :66   cara :58  
 coroa:37   coroa:32   coroa:37   coroa:34   coroa:42  

No total temos:

Exemplo 3.2

Suponhamos que temos uma caixa com 8 bolas iguais, mas com cores diferentes: 2 bolas brancas, 3 azuis, uma preta, uma vermelha e uma amarela.

A caixa com as bolas pode ser representada com ajuda do comando rep do R:

bolas <- rep(c("branca", "azul", "preta", "vermelha", "amarela"),times = c(2,3,1,1,1))
bolas

[1] "branca"   "branca"   "azul"     "azul"     "azul"     "preta"    "vermelha" "amarela" 

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

sample(bolas,4)

[1] "branca"   "preta"    "azul"     "vermelha"

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

sample(bolas,4)

[1] "branca"   "vermelha" "branca"   "preta"   

Repliquemos o experimento 10 vezes:

B <- 10
eventos <- replicate(B, sample(bolas,4))
eventos

     [,1]     [,2]       [,3]       [,4]       [,5]      [,6]      [,7]       [,8]     [,9]   
[1,] "branca" "preta"    "azul"     "amarela"  "branca"  "azul"    "vermelha" "azul"   "preta"
[2,] "azul"   "azul"     "azul"     "vermelha" "azul"    "amarela" "azul"     "branca" "azul" 
[3,] "branca" "vermelha" "preta"    "azul"     "azul"    "preta"   "azul"     "azul"   "azul" 
[4,] "preta"  "amarela"  "vermelha" "preta"    "amarela" "branca"  "amarela"  "preta"  "azul" 
     [,10]     
[1,] "amarela" 
[2,] "preta"   
[3,] "branca"  
[4,] "vermelha"

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

tab <- table(eventos)
tab

eventos
 amarela     azul   branca    preta vermelha 
       6       15        6        8        5 

Se selecionarmos 8 bolas o conjunto original é obtido:

sample(bolas,8)

[1] "branca"   "amarela"  "azul"     "azul"     "preta"    "vermelha" "branca"   "azul"    

Exemplo 3.3

Vejamos qual é a probabilidade de que, em duas escolhas, realizadas sem substituição, obtenhamos bolas azuis. Denominemos tal probabilidade por \(Pr(BB)\). É fácil mostrar que \[ Pr(BB)= \frac{3}{8}\cdot \frac{2}{7} = 0.1071429 \] Este resultado è idêntico àquele obtido no Exemplo 1.2. 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 “azul” “azul”:

experimento <- sample(bolas,2)
experimento

[1] "vermelha" "branca"  

Façamos o experimento 20 vezes:

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

 [1] FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE  TRUE FALSE
[16] FALSE FALSE FALSE FALSE FALSE

O evento BB ocorreu 2 vezes. Calculemos a proporção de número de TRUE e número de FALSE, que é dada pelo comando mean

mean(bolasazuis)

[1] 0.15

Façamos agora \(10^6\) repetições:

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

[1] 0.106618

de modo que reobtemos, aproximadamente, o resultado esperado.

Exemplo 3.4

Suponhamos que queremos saber agora qual é a probabilidade \(Pr(RWB)\) (red,white,blue) de obtermos bolas vermelha, branca e azul, respectivamente, ao retirarmos 3 bolas, sem substituição. Sabemos que \[ Pr(RWB)=\frac{1}{8}\frac{2}{7}\frac{3}{6}=\frac{1}{56}= 0.01785714285\,. \] Vejamos como obter este resultado a partir de uma simulação de Monte Carlo, fazendo \(10^6\) repetições:

B <-1000000
bolascoloridas <- replicate(B,{experimento <- sample(bolas,3) 
(experimento[1]=="vermelha" & experimento[2]=="branca" & experimento[3]=="azul")})
mean(bolascoloridas)

[1] 0.017928

que é, aproximadamente, o resultado esperado.

Exemplo 3.5

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

bolas <- rep(c("branca", "azul", "preta", "vermelha", "amarela"),times = c(2,3,1,1,1))
bolas

[1] "branca"   "branca"   "azul"     "azul"     "azul"     "preta"    "vermelha" "amarela" 

s <- sample(bolas,5)
s

[1] "branca"  "amarela" "azul"    "azul"    "branca" 

Contemos o número de bolas vermelhas na amostra:

library(stringr)
str_count(s, "azul")

[1] 0 0 1 1 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, "azul")
length(which(n==1))

[1] 2

Aplicando este comando ao nosso problema, repetindo o experimento \(10^6\) vezes, a proporção de casos nos quais temos duas ou mais bolas azuis é:

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

[1] 0.714405

Consideremos novamente o Exemplo 2.2, que agora será resolvido por meio de simulação:

Exemplo 3.6

Suponhamos que uma caixa temos 2 bolas brancas, 3 azuis, uma preta, uma vermelha e uma amarela. Qual é a probabilidade de que numa seleção de 4 bolas obtenhamos 2 bolas azuis e uma bola amarela? Solução. Representemos a caixa de bolas:

bolas <- rep(c("branca", "azul", "preta", "vermelha", "amarela"),times = c(2,3,1,1,1))
bolas

[1] "branca"   "branca"   "azul"     "azul"     "azul"     "preta"    "vermelha" "amarela" 

A proporção de casos nos quais temos duas bolas azuis e uma amarela é:

library(stringr)
B <- 1000000
bolascoloridas<-replicate(B,{experimento <- sample(bolas,4) 
n1 <- str_count(experimento, "azul")
n2 <- str_count(experimento, "amarela")
(length(which(n1==1)) ==2)   & (length(which(n2==1)) ==1)  } ) 
mean(bolascoloridas)

[1] 0.171786

que é similar ao valor 0.1714286 obtido no Exemplo 2.2.

Exemplo 3.7

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"  "blue"  "blue"  "blue"  "blue"  "blue" 
[12] "white" "white" "white" "white" "white" "white" "green" "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.1101756

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

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

[1] 0.2248067

Exemplo 3.8. Jogo de volei

Temos dois times de volei, Brasil e Argentina, que estão jogando uma série 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. Agora que as probabilidades dos eventos são diferentes, devemos especificá-las por um vetor. A opção replace = TRUE agora é necessária. A opção any de replicate significa que a condição win é verificada para todos os elementos da amostra.

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

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

Exemplo 3.9. Jogo de volei (2)

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? (i) Calcule exatamente (ii) Faça uma simulação de Monte Carlo.

Solução

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

resultados <- c(0,1)

A função list lista o conteúdo do vetor no argumento:

list(resultados)

[[1]]
[1] 0 1

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

L <- replicate(6,list(resultados))
L

[[1]]
[1] 0 1

[[2]]
[1] 0 1

[[3]]
[1] 0 1

[[4]]
[1] 0 1

[[5]]
[1] 0 1

[[6]]
[1] 0 1

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 FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[16]  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE  TRUE FALSE  TRUE
[31]  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE  TRUE FALSE
[46]  TRUE  TRUE  TRUE FALSE FALSE FALSE  TRUE FALSE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE
[61]  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.3453

de modo que obtemos um valor similar aa exato.

Exemplo 3.10 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)
p

 [1] 0.500 0.525 0.550 0.575 0.600 0.625 0.650 0.675 0.700 0.725 0.750 0.775 0.800 0.825 0.850
[16] 0.875 0.900 0.925 0.950

Definimos a função:

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

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

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

Podemos calcular uma probabilidade específica. Por exemplo,

prob_win(0.6)

[1] 0.8895

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

Consideremos um problema semelhante ao anterior, mas agora com probabilidade de que A vença B em qualquer jogo fixa em p=0.7. A variável agora é o número de jogos da série.

Solução.

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

prob_A_vence <- function(N, p=0.7){
      B <- 10000
      resultado <- replicate(B, {
        A_vence <- sample(c(1,0), N, replace = TRUE, prob = c(p, 1-p))
        sum(A_vence)>=(N+1)/2
        })
            mean(resultado)
    }

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

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

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

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

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

Referências

1. Simulations in R, Datacamp, Data Analysis and Statistical Inference https://campus.datacamp.com/courses/statistical-inference-and-data-analysis/lab-2-probability?ex=10

2. Christian Robert, George Casella. Introducing Monte Carlo Methods with R, Springer, 2010.

3. Matthias Templ.Simulation for Data Science with R Packt, 2016.

4. Vikram Dayal. Quantitative Economics with R, Springer, 2020.