Introdução

Simulação é um procedimento para aproximar probabilidades. A simulação consiste em repetir computacionalmente um grande número de vezes um particular experimento aleatório, sendo que a probabilidade do resultado de interesse é aproximada pela freqüência relativa da ocorrência desse resultado nas repetições do exeperimento. Durante a Segunda Guerra Mundial, os físicos do Laboratório de Los Alamos (desenvolveu a bomba atômica) usaram simulação para estudar sobre as características do deslocamento do nêutron através de diferentes materiais. Como este trabalho era secreto, foi usado um codinome, denominado “Monte Carlo”, em homenagem ao famoso cassino de Monte Carlo. Desde então, o uso de experimentos de simulação para entender melhor os padrões de probabilidade é denominado de Método Monte Carlo.

Funções sample e replicate:

O foco aqui será no uso das funções sample e replicate as quais simplificam bastante o processo de programação de experimentos de simulação:

Verifique atentamente a documentação dessas funções, com seus exemplos de aplicação.

Exemplos:

Jogo com moeda

Pedro e carlos jogam repetidamente uma moeda justa. Se a moeda der cara, Pedro ganha \(\$1\) de Carlos; caso contrário, Carlos ganha \(\$1\) de Pedro.

Pedro inicia com \(\$0\). O que acontece com Pedro se o jogo é disputado 50 vezes?

Em cada jogada Pedro ganha ou perde \(\$1\) com igual probabilidade. O conjunto de resultados é S <- c(1, -1)

  • Simulação de uma jogada:
S <- c(1, -1)
sample(x = S, size = 1, replace = T)
## [1] 1
  • Simulação de 50 jogadas:
(lancamentos <- sample(x = S, size = 50, replace = T))
##  [1] -1  1 -1  1 -1  1  1  1  1 -1  1  1  1 -1 -1  1 -1  1  1 -1  1 -1  1 -1  1
## [26]  1  1 -1 -1  1  1 -1  1 -1 -1  1  1  1 -1 -1 -1 -1  1  1  1  1  1 -1  1 -1
  • Ganho de Pedro em 50 jogadas:
(ganho.P <- sum(lancamentos))
## [1] 8
  • Ganho médio de Pedro por jogada, em 50 lançamentos da moeda:
(media.P <- mean(lancamentos))
## [1] 0.16
  • Histograma de ganho médio por jogada de Pedro, em 1.000 repetições de 50 jogadas:
replicas <- replicate(1000, mean(sample(S, size = 50, replace = TRUE)))
head(replicas)
## [1] -0.08  0.04 -0.12  0.08  0.08 -0.12
dim(replicas)
## NULL
hist(replicas, freq = F, ylab = "Densidade", main = "")

  • Função para construir histograma de 1.000 repetições de 50 jogadas
funcao <- function(x){
  hist(replicate(1000, mean(sample(x = S, size = 50, replace = TRUE))), 
       freq = F, ylab = "Densidade", main = "", xlab = paste("Replica #", x))
}
  • 4 histogramas de 1.000 repetições de 50 jogadas - loop com o comando lapply
par(mfrow = c(2, 2))
invisible(lapply(1:4, FUN = funcao))

Qual é a probabilidade de que Pedro termine as 50 jogadas com a mesma quantidade de dinheiro que tinha ao iniciar o jogo (\(\$0\))?
  • Ganho acumulado de Pedro nas particulares 50 jogadas armazendas em lancamentos:
(ganho.acum <- cumsum(lancamentos))
##  [1] -1  0 -1  0 -1  0  1  2  3  2  3  4  5  4  3  4  3  4  5  4  5  4  5  4  5
## [26]  6  7  6  5  6  7  6  7  6  5  6  7  8  7  6  5  4  5  6  7  8  9  8  9  8
plot(ganho.acum, type = "l", xlab = "Jogada", ylab = "Ganho acumulado ($)")
abline(h = 0)

  • Quatro trajetórias de ganho acumulado de 50 jogadas (4 particulares 50 jogadas) - loop com o comando for
par(mfrow = c(2, 2))
for(j in 1:4){
  jogadas <- sample(x = S, size = 50, replace = TRUE)
  plot(cumsum(jogadas), type = "l", ylim = c(-15, 15), xlab = "Jogada", 
       ylab = "Ganho acumulado ($)", main = paste("Trajetória #", j))
  abline(h = 0, lwd = 2)
}

  • Para obtermos aproximações de probabilidades relacionadas com ganho acumulados em 50 jogadas, devemos repetir o jogo uma grande quantidade de vezes. Para facilitar a contrução do gráfico em um loop, criaremos uma função para plotar o ganho acumulado de cada particular sequência de 50 jogadas.
funcao.acum <- function(x){
  lines(cumsum(sample(x = S, size = 50, replace = TRUE)), type = "l")
}
  • Gráfico de 1.000 trajetórias de ganho acumulado:
plot(1:50, ylim = c(-30, 30), type = "n", xlab = "Jogada", 
     ylab = "Ganho acumulado ($)")
abline(h = 0, col = "red", lwd = 2)
invisible(lapply(1:1000, funcao.acum))

Com essas simulações conseguimos aproximar várias probabilidades relacionadas com o jogo. Ex.: ganho máximo, quebra, ou ruína do jogador (quando estabelecermos um capital inicial para o jogador), quantidade de vezes que um jogador específico está à frente, etc. Todas essas probabilidades são difícieis de se calcular analiticamente. As simulações permitem uma estimação (aproximação) de seus valores.

A probabilidade solicitada é simples de ser determinada, pois é a probabilidade de Pedro terminar o jogo com \(\$0\) ou seja de ganhar (ou perder!) exatamente 25 jogadas. A variável aleatória F, associada à quantidade de caras é uma binomial com parâmentros \(n = 50\) jogadas e probabilidade de sucesso \(p = 0,5\), pois a moeda é justa. A expressão dessa probabilidade é:

\[\begin{align*} \operatorname{P}(F=0) = \binom{50}{25} 0,5^{25}0,5^{25}= \binom{50}{25} 0,5^{50} = 0,1122752 \end{align*}\]

ou, calculado diretamente pelo R por meio do comando dbinom (verifique a documentação desse comando):

dbinom(25, 50, 0.5)
## [1] 0.1122752

E para aproximarmos essa probabilidade por meio de simulação?

  • Função para simular a valor ganha (ou perdida, se a quantidade for negativa) em um jogo (n lançamentos de uma moeda honesta):
P.C <- function(n = 50){
  jogadas <- sample(x = S, size = n, replace = TRUE)
  sum(jogadas)
}

Considerado \(n= 50\) como default.

  • Simulação de um jogo com 50 (default) lançamentos de uma moeda honesta:
P.C()
## [1] -2
  • E o que acontece em termos da variável aleatória \(F\), associada ao valor ganho (ou perdido) em jogo com 50 lançamentos de uma moeda honesta:
  # Simulação com 1.000 repetições de jogo com 50 lançamentos
F <- replicate(1000, P.C())
  # Distribuição de frequência (probabilidade aproximada) dos ganhos
table(F)
## F
## -22 -20 -18 -16 -14 -12 -10  -8  -6  -4  -2   0   2   4   6   8  10  12  14  16 
##   2   1   9   8  20  22  37  72  84  97 106 102 112  97  82  60  35  25  20   4 
##  18  20  24 
##   2   2   1
  # Distribuição de frequência relativa (probabilidade aproximada) dos ganhos
prop.table(table(F))
## F
##   -22   -20   -18   -16   -14   -12   -10    -8    -6    -4    -2     0     2 
## 0.002 0.001 0.009 0.008 0.020 0.022 0.037 0.072 0.084 0.097 0.106 0.102 0.112 
##     4     6     8    10    12    14    16    18    20    24 
## 0.097 0.082 0.060 0.035 0.025 0.020 0.004 0.002 0.002 0.001
  # Gráfico da distribuição das frequências (probabilidade aproximada) de F 
plot(prop.table(table(F)), ylab = "Frequencia relativa", xlab = "Ganho final")

  • Qual a frequência relativa obtida para \(F = 0\)?
  # valor exato de P(F=0)
dbinom(25, size = 50, prob = 0.5)
## [1] 0.1122752
  # valor aproximado de P(F=0)
prop.table(table(F))[11]
##    -2 
## 0.106
  # Erro do valor aproximado
prop.table(table(F))[11] - dbinom(25, size = 50, prob = 0.5)
##           -2 
## -0.006275173

Aumentar a quantidade de réplicas melhora a aproximação? Tente.

Vamos incrementar a função, acrescentando o cálculo aproximado das vezes em que Pedro está na frente do jogo (soma acumulada positiva) e o ganho máximo durante o jogo.

P.C <- function(n = 50){
  jogadas <- sample(x = S, size = n, replace = TRUE)
  acumulado <- cumsum(jogadas)
  c(F = sum(jogadas), L = sum(acumulado > 0), M = max(acumulado))
}
  • Simulação de um jogo
P.C()
##  F  L  M 
##  2 18  3
  • Simulação de 1.000 jogos (1.000 réplicas)
S <- replicate(1000, P.C())
dim(S)
## [1]    3 1000
Qual é o número provável de lançamentos em que Pedro estará ganhando?
  # qte. vezes na frente
vezes.na.frente <- S["L", ]
t(table(vezes.na.frente))
##       vezes.na.frente
##          0   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17
##   [1,] 136  54  30  27  19  23  13  19   7  23  13  16   8  10  13   9   8  19
##       vezes.na.frente
##         18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35
##   [1,]  12  16  15  16  13  11  13  11  16  14  15  14  12  11   9  17  14  13
##       vezes.na.frente
##         36  37  38  39  40  41  42  43  44  45  46  47  48  49  50
##   [1,]  14  10  15  12   8  11  18  15  20  25  23  22  29  40  49
plot(prop.table(table(vezes.na.frente)))

Percebe-se que o jogo é justo, pois, os valores extremos são mais prováveis

Qual será o valor do ganho máximo de Pedro durante o jogo?
maximo.ganho <- S["M", ]
prop.table(table(maximo.ganho))
## maximo.ganho
##    -1     0     1     2     3     4     5     6     7     8     9    10    11 
## 0.068 0.068 0.109 0.084 0.091 0.099 0.071 0.086 0.065 0.054 0.037 0.047 0.026 
##    12    13    14    15    16    17    18    19    20    24 
## 0.032 0.011 0.025 0.005 0.007 0.002 0.008 0.003 0.001 0.001
plot(prop.table(table(maximo.ganho)))

  • Probabilidade de o ganho máximo ser maior ou igual a 10
sum(maximo.ganho >= 10)/1000
## [1] 0.168

Problema dos chapéus

Antigamente os chapéus eram muito usados e quando as pessoas iam a um restaurante deixavam seus chapéus na chapelaria, que costumava ficar próxima à porta de entrada. Suponha que \(n\) chapéus tenham sido deixados na chapelaria e que, na saída, por um erro de controle, os chapéus sejam devovildos ao acaso. O número de chapéus devolvidos aos verdadeiros donos é uma quantidade aleatória e temos interesse me estudá-la.

Vamos considerar \(n = 10\) pessoas guardando seus chapéus na entrada desse restaurante. Ele serão numerados de \(1\) a \(10\), armazendaos na variável chapeus:

n <- 10
chapeus <- 1:n

A simulação de uma devolução feita descontroladamente e ao acaso pode ser efetuada por meio do comando sample(), como, por exemplo, está apresentado abaixo. Perceba que o vetor devolucao é uma permutação do vetor chapeus:

devolucao <- sample(x = chapeus, size = n, replace = TRUE)

Verifique na documentação da função que uma permutação aleatória pode ser obtida também por meio do comando sample(chapeus).

Nossa questão de interesse é saber quantos chapéus foram devolvidos corretamente a seus verdadeiros donos. Isso é determinado facilmente pelo comando abaixo:

chapeus == devolucao
##  [1] FALSE  TRUE FALSE FALSE FALSE  TRUE  TRUE FALSE FALSE FALSE

No exemplo, o subconjunto referente à devolução correta é: {2, 6, 7}.

Por sua vez, a quantidade de chapéus devolvidos corretamente é obtida pelo comando:

n.correto <- sum(chapeus == devolucao)
n.correto
## [1] 3

Vamos construir a função mistura.chapeus, que agregue os passos acima, para utilização em um experimento mais amplo com a finalidade de calcular aproximações de probabilidades de interess. A saída dessa função é a quantidade de chapéus devolvidos corretamente a seus donos. Seu código está apresentado abaixo.

mistura.chapeus <- function(n){
   
    # ----- Argumento da função ----------
    # n: quantidade de chapéus entregues na entrada do restaurante
  
  chapeus <- 1:n
  devolucao <- sample(n)
  sum(chapeus == devolucao)
}

Vamos testar a função. Antes, por precaução, vamos remover os objetos n, chapeus e devolucao que foram criados anteriormente para verificação dos passos da função. Tome sempre o cuidado de verificar suas funções passo-a-passo e de remover objetos que tenham o mesmo nome. Assim, você conseguirá identificar rapidamente possíveis fontes de erro.

rm(n, chapeus, devolucao)

Agora podemos testar a função mistura.chapeus:

mistura.chapeus(10)
## [1] 0
mistura.chapeus(30)
## [1] 4

Claro, que cada vez que eu utilizar essa função obterei valores diferentes de chapéus devolvidos corretamente. Assim, vamos aplicar essa função uma grande quantidade de vezes para verificar a frequência dessas quantidades aleatórias de chapéus devovidos corretamente a seus donos. Uma questão que ainda não discutimos é qual a menor quantidade que poderia ser considerada grande o suficiente para produzir aproximações que consideremos boas para o problema em questão. Por enquanto, neste problema de \(10\) chapéus entregues na portaria, vamos considerar \(1.000\) uma quantidade de repetições da simulação grande o suficiente para produzir um boa aproximação das probabilidades da quantidade de chapéus devolvidos corretamente (é uma variável aleatória).

matchs <- replicate(1000, mistura.chapeus(10))
  # tabela de frequencias absolutas
absoluta <- table(matchs)
absoluta
## matchs
##   0   1   2   3   4   5 
## 377 376 186  44  13   4
  # tabela de frequencias relativas
proporcao <- prop.table(absoluta)
proporcao
## matchs
##     0     1     2     3     4     5 
## 0.377 0.376 0.186 0.044 0.013 0.004

No exemplo, A propabilidade aproximada de nenhum chapéus devolvido corretamente é 0.377.

Essas proporções serão diferente se o replicate for repetido. Lembre-se que elas são aproximações da probabilidade exata (que não calculamos). Quanto menos os valores das simulações mudarem melhor deve ser essa nossa aproximação (como não sabemos o valor exato não conseguimos calcular esse erro de aproximação). Será que as aproximações ficarão mais “estáveis” se aumentarmos a quantidade de repetições? Será que se aumentarmos a quantidade de repetições começarão a aparecer os valores extremos (todos os chapéus devolvidos corretamente, por exemplo)?

Utilizando o objeto proporcao (tabela de frequências relativas de nosso experimento) conseguimos obter aproximações de várias quantidades de interesse. Se considerarmos que esse procedimento de devolução seja estendido indefinidamente, a média (valor esperado) de chapéus devolvidos corretamente será:

mean(matchs)
## [1] 0.952

Repita o procedimento, e aproxime o valor esperado de chapéus devolvidos corretamente caso \(30\) chapéus tenham sido entregues na entrada. O que achou de interessante?

Pode-se provar que o valor esperado exato dessa variável aleatória é \(\operatorname{E}(X_n) = 1, \, \forall n\)!

E para obtermos os valores mais prováveis? Quais os valores mais comuns? Etc.

O objetivo agora é estudar mais atentamente a probabilidade de nenhum chapéu ser devolvido corretamente, ou seja, \(\operatorname{P}(X_n = 0)\). Para isso, será construída uma função para simular esse experimento \(10.000\) vezes, para então determinar a proporção de experimentos em que nenhum chapéus é devolvido corretamente.

prop.nenhum <- function(n, run = 10000){
     
    # ----- Argumentos da função ----------
    # n:   quantidade de chapéus entregues na entrada do restaurante
    # run: quantidade de repetições do experimento de devolução chapéus
  
    # ----- Saída da função ----------
    # prob: proporção de vezes com nenhum chapéus devolvido corretamente
 
 corresponde <- replicate(run, mistura.chapeus(n))
 prob <- sum(corresponde == 0)/run
 prob
}

Uma aproximação (aleatória) para \(10\) chapéus entregues na entrada é 0.372, obtida com o comando prop.nenhum(10). Por outro lado, uma aproximação para \(20\) chapéus entregues na entrada é 0.366, obtida agora com o comando prop.nenhum(20). Surge aqui uma intuição interessante: será que \(\operatorname{P}(X_n = 0)\) é a mesma para qualquer quantidade de chapéus entregues na entrada? Seria bastante similar, ao fato de que o valor esperado dessa mesma variável aleatória também ser constante (e igual a 1) para qualuqer valor de \(n\). Para explorar um pouco mais essa nossa intuição (estou esperando que seja sua intuição também), vamos visualizar graficamente as aproximações para \(\operatorname{P}(X_n = 0)\), para vários valores de \(n\) (\(2 \leq n \leq 20\))

probabilidades <- sapply(2:20, prop.nenhum)
plot(x = 2:20, y = probabilidades, ylab = "Prob(zero acerto)",
     xlab = "Quantidade de chapéus entregues na entrada")
abline(h = 0.37)

Aparentemente, a probabilidade de nenhum acerto nas devoluções parece estabilizar próxima do valor \(0,37\) para quantidades maiores de chapéus entregues na entrada.

A intuição é bastante importante na solução de problemas de probabilidade, mas deve sempre ser verificada (provada). Mãos à obra então, pessoal, provem que \(\operatorname{E}(X_n) = 1, \, \forall n\) (isso é fato) e verifiquem se \(\lim_{n \to \infty}\operatorname{P}(X_n = 0)\) converge para algum valor.

Problema do colecionador

Considere um álbum de figurinhas de futebol, uma para cada jogador de futebol importante no contexto esportivo (incluindo os jogadores da Atlética!). As figurinhas são obtidas em pacotes com 10 unidades. O objetivo do colecionador é completar o álbum (nenhuma figurinha faltando). Ao comparar pacotes com 10 unidades muitas figurinhas vêm repetidas e para completar o álbum pode ser necessário comprar muitos pacotes de figurinhas. Note que embora a quantidade de figurinhas sejas fixa (determinística) a quantidade de pacotes é variável (aleatória).

Vamos considerar que nosso álbum seja composto pelas figurinhas dos IMORTAIS jogadores da Copa de 1970, ou seja, o vetor jogadores é composto pelo nome desses fabulosos atletas:

imortais <- c("Felix", "Carlos Alberto", "Brito", "Piazza", "Everaldo", "Clodoaldo", "Gerson", "Jairzinho", "Tostão", "Pelé", "Rivelino")

Vamos simular a compra de 20 figurinhas desse maravilhoso álbum, considerando que são iguais as probabilidades de obtenção de aualquer figurinha:

compra <- sample(x = imortais, size = 20, replace = TRUE)
compra
##  [1] "Gerson"         "Felix"          "Rivelino"       "Rivelino"      
##  [5] "Everaldo"       "Clodoaldo"      "Gerson"         "Tostão"        
##  [9] "Tostão"         "Carlos Alberto" "Felix"          "Brito"         
## [13] "Tostão"         "Gerson"         "Tostão"         "Pelé"          
## [17] "Rivelino"       "Jairzinho"      "Felix"          "Felix"

Não foi comprado o conjunto completo. Quantas duplicatas foram obtidas? Quantas figurinhas faltam para completar esse álbum?

unique(compra)
##  [1] "Gerson"         "Felix"          "Rivelino"       "Everaldo"      
##  [5] "Clodoaldo"      "Tostão"         "Carlos Alberto" "Brito"         
##  [9] "Pelé"           "Jairzinho"
length(unique(compra))
## [1] 10
duplicatas <- 20 - length(unique(compra))
faltantes <- length(imortais) - length(unique(compra))

Percebe-se que há 10 duplicatas e faltam 1 figurinhas para completar o álbum.

Claro que esse problema torna-se mais complexo quando o álbum possui muitas figurinhas. Para estruturar uma simulação da compra das figurinhas vamos construir a seguinte função:

colecao <- function(n, m){
  
    # ----- Argumentos da função ----------
    # n: quantidade de figurinhas do álbum
    # m: quantidade de figurinhas compradas
  
  simula.compra <- sample(x = 1:n, size = m, replace = TRUE)
  ifelse(length(unique(simula.compra)) == n, "Sim", "Não")
}

Considere que forma adquiridas \(m = 30000\) figurinhas para completar um álbum composto por \(n = 586\) figurinhas. Assim a simulação de uma particular compra de \(3.000\) figurinhas será dada por:

colecao(n = 586, m = 3000)
## [1] "Não"

O resultado acima ocorreu para uma particular compra de \(m = 3.000\) figurinhas. Como completar ou não o álbum com esse nível de compra é aleatório, desejamos a probabilidade de completar o álbum com essa compra. Vamos simular 100 diferentes compras de \(m = 3.000\) figurinhas para encontrar um valor aproximado dessa probabilidade:

  # simulação da compra de 3.000 figurinhas
compras <- replicate(100, colecao(n = 586, m = 3000))
table(compras)
## compras
## Não Sim 
##  97   3

Ou seja, a probilidade de completar o álbum é de aproximadamente 0.03, mesmo adquirindo 5.12 vezes mais figurinhas que o total do álbum.

Suponha agora que o colecionador está indeciso assegurar-se que o álbum será completado através da compra de um número arbitrariamente grande de figurinhas. A decisão terá de ser econômica para definir a quantidade de figurinhas que deverão ser compradas. Considere que cada figurinha custa \(\$0,05\), mas que há também a oprotunidade de comprar figurinhas individuais de um negociante especializado em figurinhas difíceis. Asim, o colecionador pretende comprar as figurinhas faltantes desse negociante especializado. Vamos criar as variáveis custo.rfig e custo.nfig para os custos unitários, respectivamente, da figurinha obtida ao acaso (figurinhas compradas em envelopes fechado) e daquela adquirida diretamente com o negociante. O custo total será dado por \(custo = custo.rfig \times n.comprado + custo.nfig \times n.falta\), em que n.comprado é a quantidade de figurinhas adquiridas em envelopes e n.falta é a quantidade de figurinhas faltantes, compradas diretamente com o negociante especializado. Vamos construir uma nova função (coleca2) para agregar o custo das figurinhas no problema.

colecao2 <- function(n.comprado, n = 586){
    
    # ----- Argumentos da função ----------
    # n: quantidade de figurinhas do álbum
    # n.comprado: quantidade de figurinhas compradas

  custo.rfig <- 0.05
  custo.nfig <- 0.25
  simula.figs <- sample(1:n, size = n.comprado, replace = TRUE)
  n.fig <- length(unique(simula.figs))
  n.falta <- n - n.fig
  custo = n.comprado * custo.rfig + n.falta * custo.nfig
  custo
}

O custo para completar o álbum é aleatório e depende da quantidade de figurinhas únicas que são obtidas ao comprar n.comprado figurinhas em envelopes fechados e que varia em cada compra. Por exemplo, para uma particular compra de \(3.000\) figurinhas para completar um álbum composto por \(586\) figurinhas seu custo será:

colecao2(n.comprado = 3000, n = 586)
## [1] 151

O padrão da variável aleatória associada ao custo de comprar \(800\) figurinhas em envelopes fechados, para completar o álbum com um total de \(586\) figurinhas pode ser estimado por:

 custos <- replicate(500, colecao2(800))
summary(custos)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   71.00   76.25   77.38   77.45   78.75   83.50

Uma leitura direta permite observar que há \(50\%\) de probabilidade de o custo estar entre 76.25 e 78.75.

O razoável é basear-se no custo esperado da compra de n.comprado figurinhas ue será comprada em envelopes fechados para decidir-se pela melhor estratégia para completar esse álbum. Vamos estimar esse valor esperado através da simulação run compras. A função para a obtenção do valor aproximado do custo médio de minha compra será dada por:

custo.esperado <- function(n.comprado, run = 100){
  mean(replicate(run, colecao2(n.comprado)))
}

Para decidir-se temos que responder à pergunta: qual a quantidade de figurinhas em envelope fechado devemos adquirir para minimizar o custo esperado? Vamos utilizar o código abaixo para construir um gráfico do custo esperado para vários valores de n.comprado para obtermos um valor aproximado desde custo esperado mínimo (e, claro da quantidade de figurinhas correspondente).

N <- 500:1500
ECUSTO <- sapply(N, custo.esperado)
plot(ECUSTO ~ N, xlab = "Figurinhas compradas em envelopes fechados", 
     ylab = "Custo esperado, em pilas")
grid(col = "black")

O ideal seria enxergarmos uma linha suave por entre os pontos. Poderia ser a média, mas optou-se pelo uso de um suavizador para permitir mais liberdade no desenho dessa curva. Repete-se o código anterior usando a função loess() que é um suavizador muito utilizado em exploração de relações com dados empíricos.

N <- 500:1500
ECUSTO <- sapply(N, custo.esperado)
plot(ECUSTO ~ N, xlab = "Figurinhas compradas em envelopes fechados", 
     ylab = "Custo esperado, em pilas")
grid(col = "black")
suave <- loess(ECUSTO ~N, span = 0.2)
lines(N, predict(suave), lwd = 2, col = "white")

Com a linha suavizada, fica mais fácil observar que a quantidade de figurinhas em envelopes fechados que minimiza o custo esperado é aproximadamente \(950\), implicando num custo (também aproximado) de \(\$ 76,50\).

Você poderá encontrar mais informações sobre a função loess() na url: Loess regression with loess.

Recomendação:

Simulação é uma poderosa ferramenta para solução de problemas de probabilidade. Aproveite esse estudo dirigido e pesquise bastante além do que está sugerido. A habilidade em modelar probabilisticamente problemas, além de reforçar o entendimento detalhado da situação, permite o cálculo aproximado de quantidades (probabilidade, valor esperado, variância, etc.) que seriam bastante difíceis de serem obtidas analiticamente. Entender modelos é uma importante habilidade para os profissionais de estatística. Você irá adquirir essa habilidade treinando bastante. Não hesitem em me procurar caso tenham alguma dúvida sobre o conteúdo da disciplina.

Divirtam-se!!!.

Atividade Suplementar:

Instruções para entrega da atividade suplementar

  1. Essa atividade suplementar é opcional. Ela poderá ser considerada em substituição à qualquer atividade que não tenha sido entregue. Poderá também melhorar a nota final de acordo com sua ponderação.

  2. Esta atividade suplementar deve ser entregue no Moodle até as 23h55 do dia 13/09/20.

  3. O trabalho deve ser elaborado em formato de relatório R Markdown, extensão .html, .doc ou .pdf. Fazer também o upload do arquivo .Rmd.

  1. É muito importante que você pesquise, estude e entenda cada comando novo que for apresentado nesse estudo dirigido.

Questão 1 - Apostando na roleta

Suponha que alguém jogue roleta uma quatidade muito grande de vezes em um cassino. (Considere essa situação como sendo infinitas vezes). Se em uma única jogada, um jogador aposta $5 no “vermelho”; ele ganhará $ 5 com probabilidade 18/38 e perderá $ 5 com probabilidade 20/38, caso saia ou não saia o “vermelho”. Se esse jogo (considerada a mesma aposta) for repetido 20 vezes, então os resultados (\(+5\), quando o apostador vencer e \(-5\), quando perder) podem ser vistos como uma amostra de tamanho 20 do vetor c(5, -5), selecionada com reposição, com probabilidades dadas pelo vetor c(18/38, 20/38), respectivamente. Esses ganhos podem ser simulados usando a função sample(), com o argumento prob que estabelece as probabilidades de seleção. O comando para simular os resultados de 20 apostas consecutivas nessa roleta é dado pela expressão abaixo:

sample(x = c(5, -5), size = 20, replace = TRUE, prob = c(18/38, 20/38))
##  [1] -5 -5  5 -5 -5 -5  5 -5 -5 -5 -5 -5  5  5  5 -5  5 -5 -5  5

Pedem-se:

  1. Construa uma função para calcular a soma dos resultados (ganhos ou perdas) de 20 apostas na roleta. Use a função ´replicate()` para repetir 100 vezes esta simulação de uma sequência de 20 apostas nessa roleta. Encontre a probabilidade aproximada de que o valor final das 20 apostas seja positivo (ganho total).

  2. A quantidade de apostas vencedoras é uma variável aleatória binomial com 20 tentativas e probabilidade de sucesso igual a 18/38. Usando a função dbinom(), encontre a probabilidade exata de que o total de ganhos seja positivo e verifique quão próxima essa probabilidade exata está do valor aproximado obtido em (a).

  3. Suponha que você acompanhe seu ganho acumulado durante o jogo e registre o número de jogadas consecutivas \(P\) em que seu ganho acumulado é positivo. Se os ganho individuais são armazenados no vetor ganhos, os ganhos acumulados são obtidos por meio do comando cumsum (ganhos) e o valor de \(P\) é obtido pelo comando sum(cumsum (ganhos)> 0). Ajuste a função construída em (a) para calcular também o valor de \(P\). Simule esse processo 500 vezes e construa uma tabela de frequência dos resultados. Faça um gráfico dos resultados e discuta quais os valores de P mais prováveis de ocorrer.

Questão 2 - Qual a quantidade de alunos excluídos?

Mosteller (2010) mencionou a seguinte questão de probabilidade cuja solução pode ser obtida por meio um experimento de simulação. Suponha que haja uma turma de 10 alunos e cada aluno escolha aleatoria e independentemente dois outros membros da turma. Quantos alunos não serão escolhidos por nenhum outro aluno?

A seguir um resumo das sugestões da abordagem de Mosteller para simular esse experimento:

 estudantes <- 1:10
m <- matrix(0, nrow = 10, ncol = 10)
s <- sample(estudantes [-j], size = 2, replace = FALSE)
m[j, s] <- c (1, 1)
\[\begin{bmatrix} 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 1.00 & 0.00 & 1.00 & 0.00 \\ 0.00 & 0.00 & 1.00 & 1.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 \\ 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 1.00 & 1.00 \\ 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 1.00 & 0.00 & 0.00 & 0.00 & 1.00 \\ 0.00 & 0.00 & 0.00 & 1.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 1.00 \\ 1.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 1.00 \\ 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 1.00 & 1.00 \\ 1.00 & 1.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 \\ 0.00 & 1.00 & 0.00 & 0.00 & 0.00 & 0.00 & 1.00 & 0.00 & 0.00 & 0.00 \\ 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 1.00 & 1.00 & 0.00 \end{bmatrix}\]

Observe que a soma por coluna da matriz apresenta quais alunos não foram escolhidos por nenhum colega. (Neste exemplo, um aluno não foi escolhido.)

  1. Contrua uma função Mosteller para executar uma única simulação deste experimento.

  2. Use a função replicate() e repita essa simulação 100 vezes.

  3. Use a função table() para tabular o número de alunos não selecionados nessas 100 simulações.

  4. Qual é o número mais provável de alunos não escolhidos? Qual é a probabilidade aproximada desse resultado?

Referências

Albert, Jim, and Maria Rizzo. 2012. R by Example. Springer Science & Business Media.

Mosteller, Frederick. 2010. The Pleasures of Statistics: The Autobiography of Frederick Mosteller. Springer Science & Business Media.

R Core Team. 2015. “R: A Language and Environment for Statistical Computing.” Journal Article. http://www.R-project.org.