Sumário geral: clique aqui

Introdução

Esta nota técnica tem o objetivo didático de explorar o problema clássico do colecionador (figurinhas, bonequinhos, etc.), pela via amostral, como forma de introduzir o leitor em alguns recursos de amostragem e simulação disponíveis no software R.

1 - Descrição do problema

Suponha que uma empresa está oferecendo uma promoção colocando 1 bonequinho de superheroi nas caixas de cereais que vende. São \(k\) bonequinhos diferentes, distribuidos aleatoriamente nas caixas de forma que sempre a probabilidade de um dado bonequinho em uma caixa é \(1/k\).

Deseja-se saber:

2 - Solução teórica para situações elementares

Algumas situações dentro desse problema podem ser resolvidas com facilidade. Se comprarmos exatamente \(k\) caixas, ou seja, \(m=k\), por exemplo, temos que

É fácil entender porque isso é verdade: o número de resultados possíveis (ordem importando) nas \(k\) caixas adquiridas será \(k^k\), em que cada resultado tem a mesma chance de ocorrência. Por outro lado os resultados envolvendo obtenção dos \(k\) bonequinhos diferentes, ou seja receber os bonequinhos \(1, 2,\ldots,k\), imaginando que a ordem importa, será \(k!\).

Se \(k=5\) temos que:

## Probabilidade de 5 bonequinhos diferentes em 5 caixas
factorial(5)/5^5
## [1] 0.0384

Um gráfico que mostra a probabilidade de obtermos \(k\) bonequinhos diferentes em \(k\) caixas compradas, para \(k\) variando de \(1\) a \(20\) é apresentado a seguir:

k=1:20
Pkk=factorial(k)/k^k
plot(k,Pkk,main="Probabilidade de k bonequinhos diferentes em k caixas")

Outras inferências, com \(m>k\) podem ser mais complicadas de tratar do ponto de vista teórico (não sendo objetivo deste material tratar dessa situação).

Soluções envolvendo a esperança matemática do numero esperado de caixas necessárias para se obter os bonequinhos faltantes, dado que já se possui um certo número de bonequinhos é fácil de ser obtida, não sendo objeto desta nota técnica, que visa explorar o problema por outro caminho.

Uma solução envolvendo amostragem e simulação será desenvolvida a seguir, para ilustra conceitos de amostragem e aplicações da lei dos grandes números na estimativa de probabilidades.

3 - Solução por amostragem e simulação

A função sample(x,size,replace=TRUE) do R é usada para produzir uma amostra com tamanho definido por size, dos elementos do vetor definido por x. Se a opção replace é definida com valor FALSE (o padrão caso nada seja definido), a amostra será realizada sem reposição, para obter amostras com reposição devemos definir replace=TRUE.

Por exemplo:

## Obtendo uma amostra tamanho 10 de um vetor os números de 1 a 5
x<-c(1,2,3,4,5)  ## ou x<-1:5
sample(x,size=10,replace=TRUE)
##  [1] 3 3 4 1 2 1 1 3 4 5

Se repetir esses comandos no seu computador os valores amostrados não serão necessariamente os mesmos, dado que a amostra é aleatória e independente desta que foi realizada. Para obter o mesmo resultado apresentado neste tutorial, use preliminarmente o comando set.seed(100) antes dos comandos especificados, como é ilustrado a seguir:

## Obtendo uma amostra tamanho 5 de um vetor com 5 elementos: números de 1 a 5
set.seed(100) ## para inicializar o sorteio
x<-1:5
amostra1<-sample(x,size=10,replace=TRUE)
amostra1
##  [1] 2 2 3 1 3 3 5 2 3 1

No contexto do problema estudado, se temos \(k=5\) bonequinhos e obtemos \(m=10\) caixas de cereal, essa amostra obtida indicaria que não teríamos completado nossa coleção com os bonequinhos recebidos: falta o 4, teríamos os bonequinhos 1, 2, 3 repetidos.

Para estimar a probabilidade de tirarmos os \(k\) bonequinhos diferentes em m caixas (que teoricamente sabemos que é 0,0384, quando \(m=k=5\)) através do processo de simulação poderíamos, conceitualmente, considerar o seguinte procedimento:

O problema agora é operacionalizarmos uma forma de obter uma amostra tamanho \(n\) suficientemente grande de \(X\), e calcular \(\bar X_n\) como uma estimativa de \(\text{Pr}(C_{m,k})\).

Já sabemos como obter uma amostra de \(m\) caixas, observando os bonequinhos que foram obtidos. Precisamos de uma forma de converter esse resultado, no resultado \(X=1\) ou \(X=0\) em termos da variável Bernoulli, relacionadas aos eventos \(C_{m,k}\) e \(\bar C_{m,k}\). Isso pode ser feito facilmente com apoio da seguinte função:

## Função que retorna valor 1 conseguimos todos os números entre 1 e k
## em um vetor representado pelo argumento ``amostra''
testerep <-function(amostra,k){
  x<-intersect(amostra,1:k) ## interseção do elementos dos vetores
  if (length(x)==k) return(1) else return(0)
  ## obs: se o comprimento da interseção for k sabemos que todos os k
  ## bonequinhos foram obtidos.
}

Para observarmos a função em ação, vamos testá-la com algumas possíveis amostras e ver o resultado obtido.

## Testando a função
testerep(c(1,2,3,4,5),k=5) ## todos os bonequinhos foram obtidos
## [1] 1
testerep(c(2,3,3,4,1,9,10,3,3,2,1,5),k=10) ## sem sucesso agora 
## [1] 0

Vamos agora criar uma função para operacionalizar o processo de simulação, apresentando como resultado um vetor contendo o resultado de \(n\) simulações da nossa variável \(X\) definida acima, definindo os valores de \(X_1,X_2,\ldots,X_n\), a partir de \(n\) amostras de \(m\) caixas, em situação que há \(k\) bonequinhos diferentes.

De fato, basta que a função registre o número de situações em que \(X=1\) das \(n\) amostras, apresentando ao final a frequência de ocorrência dessa situação, correspondente a:

## calcula a frequência de sucessos (obtenção de todos 
## os bonequinhos) em n amostras de m caixas, 
## com k possíveis bonequinhos

simula<-function(n,m,k){
  sucessos<-0
  bonequinhos<-1:k
  for(i in 1:n){
    amostra<-sample(bonequinhos,m,replace=TRUE)
    x<-testerep(amostra,k) # criando resultado Bernoulli
    if (x==1) sucessos<-sucessos+1
  }
  return(sucessos/n)
}
## simulação de n=10 milhões de amostras de 5 caixas
## considerando que temos 5 bonequinhos
## Obs: pode demorar mais de 5 min para processar
##      em função do seu computador.
set.seed(100)
simula(10000000,5,5)
[1] 0.0383978

O resultado indica que \(\bar X_n=0,0383978\) é uma estimativa para \(p=\text{Pr}(C_{m,k})\), com \(k=5\) e \(m=5\). Note a boa qualidade dessa estimativa, frente ao valor teórico obtido anteriormente: \(0,0384\).

Em casos em que não conhecemos o valor ``verdadeiro’’, podemos usar o intervalo de confiança para inferências sobre esse valor, ao nível de probabilidade de interesse. Nesse caso, poderiamos inferir do próprio resultado da simulação que o verdadeiro valor estaria no intervalo \([0,03824; 0,03851]\) com 99% de probabilidade. Caso esse intervalo seja muito considerado muito largo, para uma dada aplicação, o número de simulações deverá ser aumentado. Para dividir por 2 a amplitude, teremos que aumentar 4 vezes o tamanho da amostra.

4 - Estimando probabilidades para outras situações

Podemos ter idéia mais concreta do numero \(m\) de caixas que precisaremos “comprar” para conseguirmos os \(k\) bonequinhos, com uma mínima elaboração do ferramental desenvolvido para essa finalidade nos parágrafos anteriores.

Suponha que deseja obter uma curva mostrando as probabilidades de conseguirmos \(k=5\) bonequinhos em \(m\) caixas, com \(m\) variando de \(5\) a \(20\). Certamente aumentaremos nossas chances comprando mais caixas.

## função que produz um gráfico com as 
## probabilidades de completar os k bonequinhos 
## com número de caixas variando entre m1 e m2

graficoprob<-function(n,m1,m2,k){
 prob<-0
 caixas<-m1:m2
 for(i in 1:length(caixas)){
   prob[i]=simula(n,m=caixas[i],k)
 }
 plot(caixas,prob,main="Probabilidade de completar (k=5)")
 abline(h=0.5, col="gray")
}

Podemos obter um gráfico com as estimativas das probabilidades de completar \(k=5\), com compra de caixas variando de \(5\) a \(25\), utilizando 10 mil simulações para cada número de caixas, usando:

## Testando a função (k=5, m1=5, m2=25)
set.seed(100)
graficoprob(10000,5,25,5)

O gráfico mostra que comprando 10 caixas teremos um pouco mais de 50% de probabilidade de completar os 5 bonequinhos.