Encontrei o problema Can You Pass The Cranberry Sauce? na coluna The Riddler, do site FiveThirtyEight. Uma possível tradução para ele é a seguinte:
Para celebrar o Dia de Ação de Graças, você e 19 membros da sua família estão sentados em uma mesa circular (com distanciamento social, é claro^[Foi um problema publicado em 2020, então estávamos em época de pandemia.]). Todos na mesa gostariam de uma porção de molho de cranberry, que por acaso está na sua frente no momento.
Em vez de passar o molho em um círculo, você o passa aleatoriamente para a pessoa sentada diretamente à sua esquerda ou à sua direita. Eles então fazem o mesmo, passando aleatoriamente para a pessoa à sua esquerda ou direita. Isso continua até que todos tenham, em algum momento, recebido o molho de cranberry.
Das 20 pessoas no círculo, quem tem a maior chance de ser o último a receber o molho de cranberry?
Ao conhecer este problema, imaginei que poderia ser resolvido via simulação de Monte Carlo. Minha ideia é simular 100000 vezes o cenário deste jantar e analisar os resultados dos dados produzidos.
Como o problema estava demorando muito tempo para ser resolvido no meu computador, decidi paralelizar o código para otimizar a sua execução. O pacote parallel resolve isso facilmente. Como meu computador possui 8 núcleos, abaixo estou criando um cluster com 7 núcleos, pois deixei um deles livre para outras atividades enquanto o código era executado.
A partir do for, eu começo a simulação na posição 1. A partir disso, vou sorteando a direção para a qual o molho é passado e vou atualizando estas informações. Note que, se por algum motivo a posição se tornar 0, ela automaticamente é transformada em 20. Afinal, as cadeiras estão numeradas de 1 a 20. De modo análogo, se a posição 21 for sorteada, ela automaticamente é transformada em 1. Assim, enquanto cada vetor ordem tiver menos de 20 posições distintas, a simulação é executada.
library(parallel)
clust <- makeCluster(7, type = "PSOCK")
simulacao <- list()
repl <- 100000
for (j in 1:repl){
posicao <- 1
ordem <- 1
i <- 1
while (dim(table(ordem)) < 20){
if (runif(1) < 0.5){
posicao <- posicao - 1
if (posicao == 0){
posicao <- 20
}
} else {
posicao <- posicao + 1
if (posicao == 21){
posicao <- 1
}
}
ordem[i+1] <- posicao
i <- i + 1
}
simulacao[[j]] <- ordem
}
stopCluster(clust)
Agora o objeto simulacao possui 10000 listas, onde cada uma delas corresponde a uma das simulações realizadas. Por isso, cada simulação vai ter um comprimento diferente, pois este comprimento equivale ao número de vezes que o molho foi passado entre os convidados.
A primeira análise que podemos fazer é contar o número de vezes que cada posição na mesa é acessada por último:
resultado <- unlist(lapply(simulacao, tail, 1))
resultado <- as_tibble(resultado) |>
rename(posicao = value)
resultado |>
count(posicao) |>
arrange(desc(n))
## # A tibble: 19 × 2
## posicao n
## <dbl> <int>
## 1 4 5395
## 2 7 5321
## 3 12 5314
## 4 17 5308
## 5 18 5308
## 6 6 5298
## 7 20 5272
## 8 3 5271
## 9 5 5270
## 10 11 5270
## 11 13 5265
## 12 14 5263
## 13 16 5262
## 14 19 5258
## 15 9 5247
## 16 15 5217
## 17 10 5202
## 18 8 5136
## 19 2 5123
Como é possível ver, a posição 20 foi a mais visitada, com 5353 visitas. Por outro lado, a posição 17 foi a menos visitada, com 5133 visitas, indicando uma vantagem de apenas 4,29%. Portanto, note que os valores estão relativamente próximos. Isso pode indicar que eles não são tão diferentes assim. Uma forma de medir isso é através de um gráfico, comparando as probabilidades empíricas de cada posição visitada com a probabilidade teórica de 1/19 de cada posição na mesa ficar por último:
resultado |>
count(posicao) |>
mutate(probabilidade = n/repl) |>
ggplot(aes(x = posicao, y = probabilidade)) +
geom_col() +
geom_abline(slope = 0, intercept = 1/19, linetype = "dashed", colour = "red") +
annotate(geom = "text", x = 3, y = 0.055, label = "1/19") +
scale_x_continuous(breaks = 2:20, minor_breaks = NULL) +
scale_y_continuous(breaks = seq(0, 0.06, 0.01), minor_breaks = NULL, limits = c(0, 0.06)) +
labs(x = "Lugar à mesa", y = "Probabilidade de Ficar por Último")
Outra forma seria realizar um teste de qui-quadrado para ver se as probabilidades de cada lugar à mesa são estatisticamente diferentes entre si:
chisq.test(table(resultado))
##
## Chi-squared test for given probabilities
##
## data: table(resultado)
## X-squared = 13.444, df = 18, p-value = 0.7645
Com um p-valor de 0.7645, temos uma certeza muito grande de que estas probabilidades são todas iguais entre si.
Sendo assim, é fácil ver que, surpreendentemente, todas as posições na mesa que não sejam a posição inicial (em meu caso, a posição 1) tem a mesma probabilidade de ficar por último, ou seja, 1/19.
Outra forma de resolver este problema pode ser encontrada em The ‘circular random walk’ puzzle: tidy simulation of stochastic processes in R.