Distribuição de probabilidade dos sorteios da mega sena

Introdução

A loteria de prognósticos é uma modalidade de apostas bastante difundida no mundo. Por intermédio de práticas ancestrais, anímicas e divinatórias, ao longo da história, distintos povos compartilham a crença na capacidade humana de antever a sorte.

No entanto, o interesse pelos jogos de azar contribuiu para o desenvolvimento da teoria da probabilidade, ferramenta utilizada para estimar as chances de ocorrência de um evento.

Este ramo da matemática pode ser utilizado para estudar experimentos aleatórios, isto é, fenômenos que repetidos sob as mesmas condições geralmente produzem resultados diferentes, como é o caso dos sorteios da mega sena.

Nesta modalidade de loteria, quem acertar os seis números sorteados, em um conjunto universo constituído pela sequência de 1 a 60, com incremento unitário, ganha o prêmio máximo do concurso.

Enunciado

Podemos afirmar que os números sorteados, para definir o(s) ganhador(es) da mega sena, são selecionados conforme o modelo probabilístico?

Implementação

Os dados coletados serão transformados para tratamento automatizado, submetidos a visualização gráfica e testes estatísticos para determinar a aleatoriedade da seleção.

Bibliotecas

Carregamos as seguintes bibilotecas, para adicionar funcionalidades extras:

  • tidyverse - para manipulação de dados;
  • rvest - para extração de dados da web;
  • ggplot2 - para visualização gráfica;
  • gridExtra - para ampliar as possibilidades de manipulação gráfica;
  • DT - para incluir tabelas em relatórios Rmarkdown;
  • rmdformats - temas adicionais para documentos Rmarkdown.
library(tidyverse)
library(rvest)
library(ggplot2)
library(gridExtra)
library(DT)
library(rmdformats)

Coleta de dados

A série histórica dos resultados da mega sena pode ser obtida no site http://loterias.caixa.gov.br/wps/portal/loterias/landing/megasena/ ou diretamente no link http://www1.caixa.gov.br/loterias/_arquivos/loterias/D_megase.zip:

Primeiro, criamos um objeto que armazena o caminho para acesso ao arquivo disponibilizado:

resultados <- ("http://www1.caixa.gov.br/loterias/_arquivos/loterias/D_megase.zip")

Depois de criados o diretório e arquivo temporários, baixamos o arquivo compactado:

td <- tempdir()
tf <- tempfile(tmpdir=td, fileext=".zip")
download.file(resultados, tf)

Em seguida, consultamos o caminho e o formato do arquivo descompactado:

fname <- unzip(tf, list=TRUE)$Name[1]
unzip(tf, files=fname, exdir=td, overwrite=TRUE)
fpath <- file.path(td, fname)
fpath
## [1] "/tmp/RtmpykopF5/d_mega.htm"

Tarefas de pré-processamento

Extração da tabela

Preliminarmente, utilizamos uma função que retorna todas as informações do código html e pesquisamos a quantidade de tabelas disponíveis na página :

result <- xml2::read_html(fpath)
tbls <- html_nodes(result, "table")
length((tbls))
## [1] 1

Em seguida, extraimos o conteúdo da única tabela existente:

dezenas <- as.data.frame(html_table(result, fill = FALSE))
## Warning: The `fill` argument of `html_table()` is deprecated as of rvest 1.0.0.
## An improved algorithm fills by default so it is no longer needed.

Descrição dos dados originais

Inspecionamos a estrutura básica do objeto que representa o conjunto de dados::

str(dezenas)
## 'data.frame':    2608 obs. of  21 variables:
##  $ Concurso                : int  1 2 3 3 4 5 6 7 8 9 ...
##  $ Data.Sorteio            : chr  "11/03/1996" "18/03/1996" "25/03/1996" "25/03/1996" ...
##  $ X1ª.Dezena              : int  41 9 36 36 6 1 19 56 53 55 ...
##  $ X2ª.Dezena              : int  5 39 30 30 59 19 40 38 17 43 ...
##  $ X3ª.Dezena              : int  4 37 10 10 42 46 7 21 38 56 ...
##  $ X4ª.Dezena              : int  52 49 11 11 27 6 13 20 4 54 ...
##  $ X5ª.Dezena              : int  30 43 29 29 1 16 22 3 47 8 ...
##  $ X6ª.Dezena              : int  33 41 47 47 5 2 47 5 37 60 ...
##  $ Arrecadacao_Total       : chr  "0,00" "0,00" "0,00" "0,00" ...
##  $ Ganhadores_Sena         : int  0 1 2 2 0 0 0 0 0 0 ...
##  $ Cidade                  : chr  "&nbsp" "" "" "" ...
##  $ UF                      : chr  "&nbsp" "PR" "RN" "SP" ...
##  $ Rateio_Sena             : chr  "0,00" "2.307.162,23" "391.192,51" "391.192,51" ...
##  $ Ganhadores_Quina        : int  17 65 62 62 39 98 109 100 60 17 ...
##  $ Rateio_Quina            : chr  "39.158,92" "14.424,02" "10.515,93" "10.515,93" ...
##  $ Ganhadores_Quadra       : int  2016 4488 4261 4261 3311 5399 7147 5736 5262 2175 ...
##  $ Rateio_Quadra           : chr  "330,21" "208,91" "153,01" "153,01" ...
##  $ Acumulado               : chr  "SIM" "NÃO" "NÃO" "NÃO" ...
##  $ Valor_Acumulado         : chr  "1.714.650,23" "0,00" "0,00" "0,00" ...
##  $ Estimativa_Prêmio       : chr  "0,00" "0,00" "0,00" "0,00" ...
##  $ Acumulado_Mega_da_Virada: chr  "0,00" "0,00" "0,00" "0,00" ...

Tratamento dos dados

Constatamos que a quantidade de observações diverge do número do último concurso:

tail(dezenas$Concurso, n = 1)
## [1] 2331

Recortamos um subconjunto de dados, para identificarmos alguma possível anomalia:

head(dezenas$Concurso)
## [1] 1 2 3 3 4 5

Notamos a existência de números de concursos repetidos, por isto excluimos as linhas duplicadas:

dezenas <- dezenas[!duplicated(dezenas$Concurso), ]

Inspecionamos novamente a base de dados:

str(dezenas)
## 'data.frame':    2331 obs. of  21 variables:
##  $ Concurso                : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ Data.Sorteio            : chr  "11/03/1996" "18/03/1996" "25/03/1996" "01/04/1996" ...
##  $ X1ª.Dezena              : int  41 9 36 6 1 19 56 53 55 25 ...
##  $ X2ª.Dezena              : int  5 39 30 59 19 40 38 17 43 4 ...
##  $ X3ª.Dezena              : int  4 37 10 42 46 7 21 38 56 18 ...
##  $ X4ª.Dezena              : int  52 49 11 27 6 13 20 4 54 57 ...
##  $ X5ª.Dezena              : int  30 43 29 1 16 22 3 47 8 21 ...
##  $ X6ª.Dezena              : int  33 41 47 5 2 47 5 37 60 38 ...
##  $ Arrecadacao_Total       : chr  "0,00" "0,00" "0,00" "0,00" ...
##  $ Ganhadores_Sena         : int  0 1 2 0 0 0 0 0 0 0 ...
##  $ Cidade                  : chr  "&nbsp" "" "" "&nbsp" ...
##  $ UF                      : chr  "&nbsp" "PR" "RN" "&nbsp" ...
##  $ Rateio_Sena             : chr  "0,00" "2.307.162,23" "391.192,51" "0,00" ...
##  $ Ganhadores_Quina        : int  17 65 62 39 98 109 100 60 17 251 ...
##  $ Rateio_Quina            : chr  "39.158,92" "14.424,02" "10.515,93" "15.322,24" ...
##  $ Ganhadores_Quadra       : int  2016 4488 4261 3311 5399 7147 5736 5262 2175 12590 ...
##  $ Rateio_Quadra           : chr  "330,21" "208,91" "153,01" "180,48" ...
##  $ Acumulado               : chr  "SIM" "NÃO" "NÃO" "SIM" ...
##  $ Valor_Acumulado         : chr  "1.714.650,23" "0,00" "0,00" "717.080,75" ...
##  $ Estimativa_Prêmio       : chr  "0,00" "0,00" "0,00" "0,00" ...
##  $ Acumulado_Mega_da_Virada: chr  "0,00" "0,00" "0,00" "0,00" ...

Comprovada a hipótese, reduzimos a estrutura da base de dados às dezenas sorteadas:

dz <- (dezenas[,c(3:8)])
row.names(dz) <- NULL
datatable(dz)

Análise exploratória dos dados

Espaços amostrais finitos e equiprováveis

Em uma distribuição uniforme, a probabilidade de um número isolado ser escolhido aleatoriamente é o quociente da equação P(X)=\(\frac{1}{N}\); nesse caso, P(X)=\(\frac{1}{60}\). Para eventos sem reposição, iteramos: P(Xi)=\(\frac{1}{N-(i-1)}\). Portanto, dado o espaço amostral representado pelo conjunto finito \(\Omega\) = {1, 2, 3, …, 60}, a sequência de resultados possíveis pode ser armazenada através da variável:

mega <- 1:60

O ganhador do concurso é aquele que acertar os 6 números sorteados. O reultado do sorteio é associado a uma combinação de eventos equiprováveis, supondo que todos os elementos do conjunto tenham a mesma chance de ocorrer. O número de combinações possíveis é calculado pela seguinte função:

choose(60,6)
## [1] 50063860

Tabela de frequência

Se a chance de ocorrer é a mesma para todos os números, considerando os eventos isolados sem reposição, espera-se uma distribuição uniforme dos números sorteados até o último certame.

Uma tabela de frequência ajuda a lançar luz sobre os dados, ao fornecer a contagem das observações:

names(dz) <- NULL
dz <- unlist(dz)
dzn <- as.data.frame(dz)
table(dzn)
## dzn
##   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20 
## 231 238 212 251 257 240 221 232 214 266 238 229 237 222 206 242 243 234 218 225 
##  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40 
## 203 203 254 245 215 193 250 248 240 247 219 235 256 245 242 241 250 235 220 224 
##  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60 
## 241 253 245 240 229 231 224 212 234 233 239 238 266 248 197 238 222 231 227 217

Podemos averiguar as percentagens ou proporções de cada observação em relação ao total da amostra:

prop.table(table(dzn))
## dzn
##          1          2          3          4          5          6          7 
## 0.01651652 0.01701702 0.01515802 0.01794652 0.01837552 0.01716002 0.01580152 
##          8          9         10         11         12         13         14 
## 0.01658802 0.01530102 0.01901902 0.01701702 0.01637352 0.01694552 0.01587302 
##         15         16         17         18         19         20         21 
## 0.01472901 0.01730302 0.01737452 0.01673102 0.01558702 0.01608752 0.01451451 
##         22         23         24         25         26         27         28 
## 0.01451451 0.01816102 0.01751752 0.01537252 0.01379951 0.01787502 0.01773202 
##         29         30         31         32         33         34         35 
## 0.01716002 0.01766052 0.01565852 0.01680252 0.01830402 0.01751752 0.01730302 
##         36         37         38         39         40         41         42 
## 0.01723152 0.01787502 0.01680252 0.01573002 0.01601602 0.01723152 0.01808952 
##         43         44         45         46         47         48         49 
## 0.01751752 0.01716002 0.01637352 0.01651652 0.01601602 0.01515802 0.01673102 
##         50         51         52         53         54         55         56 
## 0.01665952 0.01708852 0.01701702 0.01901902 0.01773202 0.01408551 0.01701702 
##         57         58         59         60 
## 0.01587302 0.01651652 0.01623052 0.01551552

Histograma de frequência

A contagem das observações pode ser melhor visualizada através de um histograma de frequência:

ggplot(dzn, aes(x = dz)) + 
  geom_histogram(aes(y = ..count..), bins = 300) + 
  scale_x_continuous(breaks=seq(1,60,1)) + 
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) + 
  xlab("dezenas sorteadas")

A princípio, uma rápida inspeção visual sugere o desenho de uma distribuição uniforme.

A soma dos números sorteados

Criamos uma nova coluna para armazenar a soma dos números sorteados, em cada concurso, posto que o resultado deve ser proporcional ao número de combinações possíveis. Os padrões podem ser melhor observados quanto maior o tamanho da amostra.

dezenas$Soma <- rowSums(dezenas[,c(3:8)])

O padrão da distribuição sugere que foram sorteados aleatoriamente:

h1 <- ggplot(dezenas) + 
        geom_histogram(aes(x=Soma), binwidth=10, fill="gray") + 
        xlab("soma dos números sorteados") + 
        ylab("")
h1

Experimento teórico

Podemos gerar a mesma quantidade de combinações, de acordo com o número do último concurso, e efetuar a soma dos números sorteados, para observarmos o padrão:

sortear <- function() {
  mega <- 1:60
 bolas <- sample(mega, size = 6, replace = FALSE)
      sum(bolas)
}

vezes <- last(dezenas$Concurso)

sorteios <- (replicate(vezes, sortear()))

sorteios <- as.data.frame(sorteios)

Geramos o gráfico da distribuição:

h2 <- ggplot(sorteios) + 
        geom_histogram(aes(x=sorteios), binwidth=10, fill="gray") + 
        xlab("soma dos números simulados") + 
        ylab("")
h2

Dispostos lado a lado, podemos comparar graficamente o comportamento dos dados observados nos sorteios com o comportamento dos dados do experimento teórico:

grid.arrange(h1, h2, ncol=2)

Teste Qui-Quadrado

Se ainda assim hesitamos em determinar se os números foram selecionados de acordo com o modelo não-determinístico, aplicamos o Teste Qui-Quadrado, utilizado para testar a aderência de um modelo probabilístico a um conjunto de dados observados:

soma <- dezenas$Soma
chisq.test(soma)
## 
##  Chi-squared test for given probabilities
## 
## data:  soma
## X-squared = 20575, df = 2330, p-value < 2.2e-16

Conclusão

2.2e-16 é a notação científica de 0.00000000000000022. Como o valor de p tende a zero, a hipótese de que os números sorteados na mega sena não sejam escolhidos aleatoriamente deve ser rejeitada.

Referências:

Grolemund, G. (2014). Hands-On Programming with R.
Ross, K. (2020). An Introduction to Probability and Simulation.
Templ, M. (2016). Simulation for Data Science with R.