Última atualização: 12 de maio, 2025 às 20:38.

1 Probabilidades

O objetivo desta lista é ilustrar a aproximação da probabilidade de eventos por simulação.
Em cada exercício, deve-se tentar encontrar, quando possível teoricamente, a probabilidade de eventos ou quantidades solicitadas.
Na sequencia, as probabilidades ou quantidades devem ser encontradas por simulação.

1.1 Encontre as probabilidades a seguir:

  • No lançamento de um dado:
    • sair a face 5;
    • sair a face 2 ou 5.
  • No lançamento de dois dados:
    • a soma dos valores ser 11;
    • a soma dos valores ser 7;
    • alguma face ser 5;
    • o valor esperado da soma.
  • No lançamento de três dados:
    • a soma dos valores ser 18;
    • a soma dos valores ser 14, 15 ou 16;
    • alguma face ser 5;
    • sair ao menos duas faces 6;
    • o valor esperado da média dos valores.
N <- 1E6
dado1 <- sample(1:6, N, replace = TRUE)
typeof(dado1)
# [1] "integer"
mean(dado1 == 5)   
# [1] 0.166677
1/6
# [1] 0.1666667
## ...
dado2 <- sample(1:6, N, replace = TRUE)
dado3 <- sample(1:6, N, replace = TRUE)
## ...
mean(dado1 == 5 | dado2 == 5 | dado3 == 5)
# [1] 0.421441
Ind5 <- (dado1 == 5) + (dado2 == 5) + (dado3 == 5)
mean(Ind5 > 0)
# [1] 0.421441
## ...

Comandos alternativos para obter soluções e comentários:

D2 <- matrix(sample(1:6, size = 2*N, replace = TRUE), ncol = 2)
dim(D2)
# [1] 1000000       2
S2 <- apply(D2, 1, sum) ## mais geral, mas ineficiente neste caso
S2 <- rowSums(D2)       ## recurso disponível para médias e somas
##
mean(S2 == 11)
# [1] 0.055351
mean(S2 == 7)
# [1] 0.16645
prop.table(table(S2))
# S2
#        2        3        4        5        6        7        8        9 
# 0.027783 0.055830 0.083751 0.110781 0.138596 0.166450 0.139146 0.111304 
#       10       11       12 
# 0.083342 0.055351 0.027666
rbind(prop.table(table(S2)), c(1:6, 5:1)/36)
#               2          3          4         5         6         7         8
# [1,] 0.02778300 0.05583000 0.08375100 0.1107810 0.1385960 0.1664500 0.1391460
# [2,] 0.02777778 0.05555556 0.08333333 0.1111111 0.1388889 0.1666667 0.1388889
#              9         10         11         12
# [1,] 0.1113040 0.08334200 0.05535100 0.02766600
# [2,] 0.1111111 0.08333333 0.05555556 0.02777778
## tempo de execução: apply vs operação vetorial
system.time({
    fun <- function(x) x[1] == 5 | x[2] == 5 
    mean(apply(D2, 1, fun))
})
#   usuário   sistema decorrido 
#     0.915     0.001     0.916
system.time({D2[,1] == 5 |  D2[,2] == 5})
#   usuário   sistema decorrido 
#     0.006     0.002     0.008

1.2 Uma oficina realiza, sequencialmente, dois serviços de rotina em seus atendimentos (S1 e S2). O tempo para realizar os serviços têm distribuição uniforme (contínua), sendo o primeiro serviço entre 15 e 25 minutos e o segundo entre 30 a 50 minutos.

  1. Qual a probabilidade do primeiro serviço:
  1. ser concluído em menos de 22 minutos?
  2. ser concluído entre 18 e 20 minutos?
  1. Qual a probabilidade do segundo serviço:
  1. ser concluído em menos de 35 minutos?
  2. ser concluído em mais de 42 minutos?
  3. ser concluído em menos de 42 minutos, sabendo-se que levou mais que 40 minutos?
  1. Para o tempo total dos serviços:
  1. Qual a probabilidade de ambos serviços serem concluídos em mais de uma hora?
  2. E entre 57 e 63 minutos?
  3. Se sabe-se que levou mais de uma hora, qual a probabilidade de ter sido alem de 70 minutos?
  4. Qual o tempo esperado de conclusão?
  1. Faça gráficos das distribuições de \(S1\), \(S2\) e \(S1+S2\).
## soluções exatas
## a.1
punif(22, min = 15, max = 25)
# [1] 0.7
## a.2
punif(20, min = 15, max = 25) - punif(18, min = 15, max = 25)
# [1] 0.2
## ou ...
diff(punif(c(18, 20), min = 15, max = 25))
# [1] 0.2
## b.3
diff(punif(c(40, 42), min = 30, max = 50))/punif(40, min = 30, max = 50, lower = FALSE)
# [1] 0.2
## aproximações por simulação
N <- 1E6
S1 <- runif(N, min = 15, max = 25)
S2 <- runif(N, min = 30, max = 50)
S12 <- S1 + S2
## a.1
mean(S1 < 22)   
# [1] 0.699867
## a.2
mean(S1 > 18 & S1 < 20)
# [1] 0.199544
## b.3
mean(subset(S2, S2 > 40) < 42)
# [1] 0.1997632
## c.2
mean(S12 > 57 & S12 < 63)
# [1] 0.300262

1.3 Considere o contexto do problema anterior porém com a seguinte condição: se S1 demora até 18 minutos, S2 tem tempo uniforme entre 30 e 50 minutos. Porém, se S1 demora mais que 18 minutos, S2 tem distribuição uniforme entre 40 a 60 minutos.

  1. Qual a probabilidade do primeiro serviço:
  1. ser concluído em menos de 22 minutos?
  2. ser concluído entre 18 e 20 minutos?
  1. Qual a probabilidade do segundo serviço:
  1. ser concluído em menos de 35 minutos?
  2. ser concluído em mais de 42 minutos?
  3. ser concluído em menos de 42 minutos, sabendo-se que levou mais que 40 minutos?
  1. Para o tempo total dos serviços:
  1. Qual a probabilidade de ambos serviços serem concluídos em mais de uma hora?
  2. E entre 57 e 63 minutos?
  3. Se sabe-se que levou mais de uma hora, qual a probabilidade de ter sido alem de 70 minutos?
  4. Qual o tempo esperado de conclusão?
  1. Faça gráficos das distribuições de \(S1\), \(S2\) e \(S1+S2\).
par(mfrow=c(1,3), mar = c(3,3.5,1,0.5), mgp = c(2,1,0))
hist(S1, nclass = 30)
hist(S2, nclass = 30)
hist(S12, nclass = 30)

1.4 Seja \(Y\) uma v.a. com distribuição normal \(Y \sim {\rm N}(\mu = 100, \sigma^2 = 10^2)\). Calcule de forma “exata” e aproximações por simulação:

  1. \(P(Y < 108)\),
  2. \(P(Y > 112)\),
  3. \(P(90 < Y < 120)\),
  4. o valor de \(a\) tal que \(P(Y < a) = 0.80\),
  5. o valor de \(a\) tal que \(P(Y > a) = 0.60\),
  6. \(P(Y < 102| Y < 112)\),
  7. \(P(Y > 115| Y > 105)\),
  8. \(P(Y < 105| Y > 92)\),
  9. \(P(100 < Y < 120| Y > 105)\).
N <- 1E6
simN <- rnorm(N, m = 100, sd = 10)
## a)
pnorm(108, m=100, sd = 10)
# [1] 0.7881446
mean(simN < 108)
# [1] 0.787664
## c)
diff(pnorm(c(90, 120), m=100, sd = 10))
# [1] 0.8185946
mean(simN > 90 & simN < 120)
# [1] 0.818972
## d)
qnorm(0.80, m=100, sd = 10)
# [1] 108.4162
quantile(simN, probs = 0.80)
#      80% 
# 108.4298
## i)
diff(pnorm(c(105, 120), m=100, sd = 10))/pnorm(105, m = 100, sd = 10, lower = FALSE)
# [1] 0.9262646
simNM105 <- subset(simN, simN > 105)
mean(simNM105 > 100 & simNM105 < 120)
# [1] 0.9259143

Respostas:

#       teorico    simulado
# a   0.7881446   0.7881390
# b   0.1150697   0.1154030
# c   0.8185946   0.8186810
# d 108.4162123 108.4153421
# e  97.4665290  97.4675945
# f   0.6545823   0.6539113
# g   0.2165286   0.2164068
# h   0.6085267   0.6079540
# i   0.9262646   0.9261101

1.5 Considere um jogo que consiste em várias rodadas independentes com probabilidade de \(0,35\) de ganhar uma rodada. O jogo termina quando se perde a terceira rodada. Calcule:

  1. Qual a probabilidade de conseguir 2 vitórias antes do jogo encerrar?
  2. Qual a probabilidade do jogo ter quatro rodadas?
  3. Qual a probabilidade do jogo ter mais de cinco rodadas?
  4. Qual o número esperado de vitórias?
  5. Qual o número esperado de rodadas?
  6. Repita os cálculos para ao menos três outros valores da probabilidade de ganhar uma rodada e número de derrotas que determina o encerramento do jogo.
  7. Faça um gráfico do número esperado de rodadas em função do valor da probabilidade.

Dica: procure identificar uma distribuição de probabilidades relacionada ao problema.

Solução: \[ \begin{align*} R &: \mbox{número de rodadas no jogo} \\ Y &: \mbox{número de vitórias até o encerramento do jogo} (Y = R - 3)\\ Y &\sim {\rm BN}(r = 3, p = 1-0.35 = 0.65) \\ P[Y=y] &= \binom{y+r-1}{y} p^r (1-p)^y \\ {\rm E}(Y) &= r \frac{1-p}{p} \mbox{ e } {\rm V}(Y) = r \frac{1-p}{p^2} \\ a)& P[Y = 2] = 0.20185\\ b)& P[R = 4] = P[Y = 2] = 0.28836\\ c)& P[R > 5] = P[Y > 2] = 0.23517\\ d)& {\rm E}(Y) = 1.615\\ e)& {\rm E}(R) = 4.615\\ \end{align*} \]

Soluções (teóricas) pela distribuição binomial negativa.

pwin <- 0.35
p <- 1-pwin  ## p é a probabilidade associada ao evento de parada
r <- 3
choose((2+r-1), 2) * (1-p)^2 * p^3
# [1] 0.2018494
##
(qa <- dnbinom(2, size = r, prob = p))
# [1] 0.2018494
(qb <- dnbinom((4-r), size = r, prob = p))
# [1] 0.2883562
(qc <- 1-pnbinom((5-r), size = r, prob = p))
# [1] 0.2351694
(qd <- r * (1-p)/p)
# [1] 1.615385
(qe <- r + r * (1-p)/p)
# [1] 4.615385

Aproximação por simulação de Monte Carlo.
Inicialmente define-se uma função para uma rodada do jogo.
Depois de testar a função simula-se jogadas consecutivas utilizando replicate().

game <- function(pwin, nlosses, mensagem = FALSE){
    nl <- nr <- 0
    if(mensagem) cat("Sequencia no jogo: ")
    while(nl < nlosses){
        u <- runif(1)
        perdeu <- u > pwin
        if(mensagem) cat(as.numeric(perdeu))
        if(perdeu) nl <- nl + 1
        nr <- nr + 1
    }
    if(mensagem) cat("\n")
    return(nr)
}
game(0.35, 3, mensagem = TRUE)
# Sequencia no jogo: 111
# [1] 3
game(0.70, 3, mensagem = TRUE)
# Sequencia no jogo: 1000101
# [1] 7

replicate(10, game(0.35, 3, mensagem = TRUE))
# Sequencia no jogo: 111
# Sequencia no jogo: 0111
# Sequencia no jogo: 0111
# Sequencia no jogo: 0100011
# Sequencia no jogo: 00111
# Sequencia no jogo: 0111
# Sequencia no jogo: 111
# Sequencia no jogo: 100011
# Sequencia no jogo: 01101
# Sequencia no jogo: 1011
#  [1] 3 4 4 7 5 4 3 6 5 4
replicate(10, game(0.70, 3, mensagem = TRUE))
# Sequencia no jogo: 01101
# Sequencia no jogo: 00010101
# Sequencia no jogo: 0010000100000001
# Sequencia no jogo: 0111
# Sequencia no jogo: 00001011
# Sequencia no jogo: 01001001
# Sequencia no jogo: 10000011
# Sequencia no jogo: 00000010011
# Sequencia no jogo: 0010100001
# Sequencia no jogo: 0000000010100001
#  [1]  5  8 16  4  8  8  8 11 10 16

Tudo parece estar OK.  “Joga-se” agora um grande número de vezes para responder questões com aproximações por simulação.

N <- 1E6
sim <- replicate(N, game(pwin, r))
mean((sim - r) == 2) 
# [1] 0.202281
mean(sim == 4)
# [1] 0.28871
mean(sim > 5)
# [1] 0.23543
mean(sim) - 3
# [1] 1.616698
mean(sim)
# [1] 4.616698

Cálculos por simulação podem se beneficiar de paralelização. Há vários recursos para paralelização, como por exemplo o pacote mcreplicate para paralelizar loops deste tipo. O código a seguir mostra como obter os tempos computacionais sem e com paralelização. Os tempos irão variam com processadores e número de núcleos de processamento, portanto o importante não são os valores obtidos dos tempos individualmente, mas a comparação destes.

now <- Sys.time()
sim <- replicate(N, game(pwin, r))
Sys.time() - now
##
now <- Sys.time()
sim <- mcreplicate::mc_replicate(N, game(pwin, r))
Sys.time() - now

Agora o gráfico relacionando o número esperado de rodadas com diferentes valores da probabilidade de ganhar uma rodada do jogo.

Resp.teo.f <- function(pwin, r){
    p <- 1 - pwin
    r + r * (1-p)/p
}
curve(Resp.teo.f(x, r = 3), from = 0.05, 0.95,
      xlab = "pwin", ylab = "N esperado de rodadas")
##
Resp.sim.f <- function(pwin, r = 3, N=1E5){
    mean(replicate(N, game(pwin, r)))
}
pseq <- seq(0.05, 0.95, by = 0.05)
Rsim <- sapply(pseq, Resp.sim.f, r = 3)
points(pseq, Rsim, col = 4, pch = 19)
legend("topleft", c("teórico", "por simulação"),
       col = c(1,4), pch = c(3, 19))

1.6 Considere o jogo anterior mas com um diferente critério de parada. O jogo se encerra quando se perde duas rodadas consecutivas. Calcule:

  1. Qual a probabilidade de conseguir duas vitórias antes do jogo encerrar?
  2. Qual a probabilidade do jogo ter quatro rodadas?
  3. Qual a probabilidade do jogo ter mais de cinco rodadas?
  4. Qual o número esperado de vitórias?
  5. Qual o número esperado de rodadas?
  6. Repita os cálculos para ao menos três outros valores da probabilidade de ganhar uma rodada e número de derrotas que determina o encerramento do jogo.
  7. Faça um gráfico do número esperado de rodadas em função do valor da probabilidade.

1.7 Probabilidades de Chevalier

Atribui-se o início do tratamento matemático formal da teoria de probabilidades a problemas envolvendo jogos.
Em particular é citado o problema de Chevalier de Méré [@spiegelhalter2019art, Capítulo 8] dos anos de 1650’s, em que buscava-se saber em qual das duas situações a seguir haveria a maior chance de ganhar.

  1. Jogo 1: joga-se um dado (honesto) no máximo quatro vezes e vence se obtiver pelo menos um seis.
  2. Jogo 2: joga-se dois dados (honestos) no máximo 24 vezes e vence se obtiver uma dupla de seis.

Chevalier de Méré decidiu então, seguindo os melhores princípios, jogar cada uma dos jogos inúmeras vezes e verificar com que frequencia ganhava em cada caso.
Baseado neste experimento ele verificou qual em qual dos jogos tinha melhores chances.
Ele não tinha, à época, o conhecimento para fazer cálculos exatos destas probabilidades.

Vamos explorar este problema para ilustrar como algorítmos computacionais podem simular experimentos que permitem avaliar probabilidades de interesse.
Neste caso, é possível obter o resultado com uma solução formal (analítica). Entretanto, a solução computacional é mais geral, portanto útil para tratar outros problemas que não podem ser resolvidos analiticamente.

  1. Obtenha inicialmente a solução formal (analítica) do problema calculando a probabilidade de vencer em cada jogo.
  2. Repita o experimento prático de Chevalier, porém usando uma simulação computacional para obter a probabilidade de vencer em cada jogo.
  3. Detalhe os resultados do experimento computacional mostrando como as probabilidades se comportam para números crescentes de experimentos computacionais.

Solução parcial

Sejam os eventos:
\(J1\): Obter pelo menos um 6 em quatro lançamentos de um dado.
\(J2\): Obter pelo menos um par de 6 em vinte quatro lançamentos de um par de dados.

Para \(J1\) a probabilidade é de se obter o 6 no primeiro ou segundo ou terceiro ou quarto lançamento, portanto: \[P[J1] = \frac{1}{6} + \left(\frac{5}{6}\right)\frac{1}{6} + \left(\frac{5}{6}\right)^2\frac{1}{6} + \left(\frac{5}{6}\right)^3\frac{1}{6} = 0.5177.\] Alternativamente e de forma mais direta pode-se calcular a probabilidade do evento complementar \(\overline{J1}\) (não obter nenhum 6 em quatro lançamentos e subtrair de 1. \[P[J1] = 1 - P[\overline{J1}] = 1 - \left(\frac{5}{6}\right)^4 = 0.5177.\]

De forma análoga para \(J2\) a probabilidade é de se obter o par de 6 no primeiro ou seguindo ou terceiro ou quarto até o décimo quarto lançamento, portanto lembrando que a probabilidade de um par de 6 é de 1/36 temos: \[P[J2] = \frac{1}{36} + \left(\frac{35}{36}\right)\frac{1}{36} + \left(\frac{35}{36}\right)^2\frac{1}{36} + \ldots + \left(\frac{35}{36}\right)^{23}\frac{1}{36} = 0.4914.\] Alternativamente e de forma mais direta pode-se calcular a probabilidade do evento complementar \(\overline{J1}\) (não obter nenhum 6 em quatro lançamentos e subtratir de 1. \[P[J2] = 1 - P[\overline{J2}] = 1 - \left(\frac{35}{36}\right)^{24} = 0.4914.\]

Solução analítica calculada (de diferentes formas) com comandos do \(\textsf{R}\).

Emulando uma rodada de cada jogo temos:

## Jogo 1:
(simJ1 <- sample(1:6, 4, replace = TRUE))
# [1] 3 6 6 1
6 %in% simJ1 
# [1] TRUE
## Jogo 2:
(simJ2 <- matrix(sample(1:6, 2*24, replace = TRUE), ncol = 2))
#       [,1] [,2]
#  [1,]    5    3
#  [2,]    5    1
#  [3,]    4    2
#  [4,]    1    1
#  [5,]    2    1
#  [6,]    6    1
#  [7,]    5    5
#  [8,]    5    6
#  [9,]    1    6
# [10,]    3    2
# [11,]    2    1
# [12,]    2    4
# [13,]    4    3
# [14,]    1    2
# [15,]    3    5
# [16,]    3    5
# [17,]    6    4
# [18,]    6    4
# [19,]    3    6
# [20,]    1    5
# [21,]    4    3
# [22,]    1    5
# [23,]    4    1
# [24,]    1    3
any(rowSums(simJ2) == 12)
# [1] FALSE

Agora é com você explorar simulações!

1.8 Um exemplo com IMC

O IMC (índice de massa corporal) usa altura e peso de indivíduos segundo a equação: \[ \text{IMC} = \frac{peso (kg)}{altura^2 (m)} . \] Segundo este site o valor obtido é utilizado para classificar o indivíduo em uma das seguintes categorias:

  • 18,5 ou menos: Abaixo do normal
  • Entre 18,6 e 24,9: Normal
  • Entre 25,0 e 29,9: Sobrepeso
  • Entre 30,0 e 34,9: Obesidade grau I
  • Entre 35,0 e 39,9: Obesidade grau II
  • Acima de 40,0: Obesidade grau III
  1. Suponha que uma população tem distribuição normal com valor médio de IMC igual a 27 e variância 25. Simule 100 valores desta população e agrupe nas classes de IMC tabulando quantos foram simulados em cada classe. Dica: cut().
  2. Suponha que uma população tem dstribuição gama com valor médio de IMC igual a 27 e variância 25. Simule 100 valores desta população e agrupe nas classes de IMC tabulando quantos foram simulados em cada classe. Dica: rgamma().
    OBS Os valores de média (27) e variância (25) informados correspondem a escolhas de \(shape = \frac{27^2}{25}\) e \(scale = \frac{27}{25}\) para usar em rgamma().

A Figura a seguir ilustra o problema. As probabilidade indicadas são valores teóricos das distribuições normal (usando pnorm()) e gamma (usando pgamma()). Compare com os valores obtidos por simulação.

amN <- rnorm(1000, mean=27, sd=5)
mean(amN)
# [1] 26.95239
var(amN)
# [1] 22.71763
##
cortes <- c(-Inf, 18.6, 25, 30, 35, 40, Inf)
classes <- c("Abaixo do normal","Normal","Sobrepeso","Obesidade I","Obesidade II","Obesidade III")
## diferentes opções para tabela:
table(cut(amN, breaks = cortes))
# 
# (-Inf,18.6]   (18.6,25]     (25,30]     (30,35]     (35,40]   (40, Inf] 
#          44         288         408         218          40           2
table(cut(amN, breaks = cortes, right=FALSE))
# 
# [-Inf,18.6)   [18.6,25)     [25,30)     [30,35)     [35,40)   [40, Inf) 
#          44         288         408         218          40           2
table(cut(amN, breaks = cortes, labels = classes))
# 
# Abaixo do normal           Normal        Sobrepeso      Obesidade I 
#               44              288              408              218 
#     Obesidade II    Obesidade III 
#               40                2
prop.table(table(cut(amN, breaks = cortes, labels = classes)))
# 
# Abaixo do normal           Normal        Sobrepeso      Obesidade I 
#            0.044            0.288            0.408            0.218 
#     Obesidade II    Obesidade III 
#            0.040            0.002