Última atualização: 12 de maio, 2025 às 20:38.
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.
1E6
N <- sample(1:6, N, replace = TRUE)
dado1 <-typeof(dado1)
# [1] "integer"
mean(dado1 == 5)
# [1] 0.166677
1/6
# [1] 0.1666667
## ...
sample(1:6, N, replace = TRUE)
dado2 <- sample(1:6, N, replace = TRUE)
dado3 <-## ...
mean(dado1 == 5 | dado2 == 5 | dado3 == 5)
# [1] 0.421441
(dado1 == 5) + (dado2 == 5) + (dado3 == 5)
Ind5 <-mean(Ind5 > 0)
# [1] 0.421441
## ...
Comandos alternativos para obter soluções e comentários:
matrix(sample(1:6, size = 2*N, replace = TRUE), ncol = 2)
D2 <-dim(D2)
# [1] 1000000 2
apply(D2, 1, sum) ## mais geral, mas ineficiente neste caso
S2 <- rowSums(D2) ## recurso disponível para médias e somas
S2 <-##
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({
function(x) x[1] == 5 | x[2] == 5
fun <-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
## 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
1E6
N <- runif(N, min = 15, max = 25)
S1 <- runif(N, min = 30, max = 50)
S2 <- S1 + S2
S12 <-## 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
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)
1E6
N <- rnorm(N, m = 100, sd = 10)
simN <-## 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
subset(simN, simN > 105)
simNM105 <-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
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.
0.35
pwin <- 1-pwin ## p é a probabilidade associada ao evento de parada
p <- 3
r <-choose((2+r-1), 2) * (1-p)^2 * p^3
# [1] 0.2018494
##
dnbinom(2, size = r, prob = p))
(qa <-# [1] 0.2018494
dnbinom((4-r), size = r, prob = p))
(qb <-# [1] 0.2883562
1-pnbinom((5-r), size = r, prob = p))
(qc <-# [1] 0.2351694
r * (1-p)/p)
(qd <-# [1] 1.615385
r + r * (1-p)/p)
(qe <-# [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().
function(pwin, nlosses, mensagem = FALSE){
game <- nr <- 0
nl <-if(mensagem) cat("Sequencia no jogo: ")
while(nl < nlosses){
runif(1)
u <- u > pwin
perdeu <-if(mensagem) cat(as.numeric(perdeu))
if(perdeu) nl <- nl + 1
nr + 1
nr <-
}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.
1E6
N <- replicate(N, game(pwin, r))
sim <-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.
Sys.time()
now <- replicate(N, game(pwin, r))
sim <-Sys.time() - now
##
Sys.time()
now <- mcreplicate::mc_replicate(N, game(pwin, r))
sim <-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.
function(pwin, r){
Resp.teo.f <- 1 - pwin
p <-+ r * (1-p)/p
r
}curve(Resp.teo.f(x, r = 3), from = 0.05, 0.95,
xlab = "pwin", ylab = "N esperado de rodadas")
##
function(pwin, r = 3, N=1E5){
Resp.sim.f <-mean(replicate(N, game(pwin, r)))
} seq(0.05, 0.95, by = 0.05)
pseq <- sapply(pseq, Resp.sim.f, r = 3)
Rsim <-points(pseq, Rsim, col = 4, pch = 19)
legend("topleft", c("teórico", "por simulação"),
col = c(1,4), pch = c(3, 19))
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.
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.
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:
sample(1:6, 4, replace = TRUE))
(simJ1 <-# [1] 3 6 6 1
6 %in% simJ1
# [1] TRUE
## Jogo 2:
matrix(sample(1:6, 2*24, replace = TRUE), ncol = 2))
(simJ2 <-# [,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!
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:
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.
rnorm(1000, mean=27, sd=5)
amN <-mean(amN)
# [1] 26.95239
var(amN)
# [1] 22.71763
##
c(-Inf, 18.6, 25, 30, 35, 40, Inf)
cortes <- c("Abaixo do normal","Normal","Sobrepeso","Obesidade I","Obesidade II","Obesidade III")
classes <-## 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