Motivação

Recentemente encontrei o artigo How to find pi in randomness all around you. Nele, a autora cita um paper intitulado Estimating π with a Coin, no qual o autor propõe como estimar o valor de π utilizando o lançamento de moedas honestas.

Segundo o paper, lançar uma moeda honesta um número grande de vezes é capaz de estimar o valor de π. Segundo o próprio texto,

Toss a coin until the first time the observed proportion of heads exceeds ½ (that is, until the cumulative number of heads exceeds the cumulative number of tails) and record that fraction; for instance, if the tosses proceeded as tails, heads, tails, heads, heads you would record the fraction ⅗. Now perform a second run of coin-tosses, starting afresh and again recording a fraction that exceeds ½ (starting afresh means that you ignore all earlier tosses). Repeat. As you perform more and more such trials, the average of the fractions will approach π⁄4.

Decidi implementar isso computacionalmente para verificar este fato, apesar do partigo demonstrá-lo matematicamente.

Função para Simulação

A função f_lancamento implementa cada sequência particular de lançamentos até a obtenção de uma sequência de caras cuja proporção em relação ao total de lançamentos seja estritamente maior do que 0.5. Neste caso, estou considerando a realização de uma variável aleatória binomial com \(n = 1\) e \(p = 0.5\) em que o resultado 1 representa cara e 0 representa coroa.

f_lancamento <- function(){
  
  lancamento <- rbinom(1, 1, 0.5)
  
  if (lancamento == 1){
    razao <- 1
  } else {
    lancamento <- c(lancamento, rbinom(1, 1, 0.5))
    while(mean(lancamento) <= 0.5){
      lancamento <- c(lancamento, rbinom(1, 1, 0.5))
    }
  }
  
  return(lancamento)

}

Abaixo é possível ver quatro realizações desta função:

set.seed(423809)

f_lancamento()

## [1] 1

f_lancamento()

## [1] 0 1 0 1 1

f_lancamento()

## [1] 1

f_lancamento()

##  [1] 0 0 0 1 0 1 1 1 0 1 0 0 0 1 0 1 0 1 0 1 1 1 0 0 0 1 0 0 1 0 1 1 1 1 1

Entretanto, simulações grandes deste tipo podem gerar problemas de gerenciamento de memória, já que a sequência inteira de lançamentos das moedas fica armazenada na memória do R. Sendo assim, a função f_lancamento_media já retorna os valores da proporção de caras em relação ao total de lançamentos da moeda.

f_lancamento_media <- function(){
    
  lancamento <- rbinom(1, 1, 0.5)
  
  if (lancamento == 1){
    razao <- 1
  } else {
    lancamento <- c(lancamento, rbinom(1, 1, 0.5))
    while(mean(lancamento) <= 0.5){
      lancamento <- c(lancamento, rbinom(1, 1, 0.5))
    }
  }
  
  return(mean(lancamento))
}

Simulação de Fato

Com a função a ser utilizada na simulação implementada, basta rodá-la um número grande de vezes para obter a estimativa de π. Como é algo computacionalmente intenso, opter por repeti-la apenas 500 vezes e paralelizar o código utilizando todos menos um núcleo do meu computador.

library(foreach)
library(doParallel)

# configurar cluster com o numero maximo de núcleos
# em meu computador, menos 1

cl <- makeCluster(detectCores()-1)
registerDoParallel(cl)

n <- 500

set.seed(1)

# loop paralelizado com foreach

system.time(

  serie <- foreach(j = 1:n, .combine = c) %dopar% {
    mean(f_lancamento_media())
  }

)

##    user  system elapsed 
##   0.060   0.016 291.392

mean(serie) * 4

## [1] 3.151777

stopCluster(cl)

Como é possível ver, a estimativa de π ficou igual a 3.1518, uma diferença de -0.0102 em relação ao valor nomimal de π. Não é uma má aproximação, além de ter sido relativamente rápida.

Entretanto, se aumentarmos o número de replicações, o tempo de execução aumenta de maneira não-linear. Afinal, ao fazer isso, aumentamos a chance de termos sequências mais longas, que demoram mais a ter uma proporção maior de caras em relação a coroas.

No fim, esta forma de cálculo de π acaba sendo mais uma curiosidade do que algo a ser utilizado de forma prática. Entretanto, é interessante que apenas em 2025 alguém tenha conseguido pensar nela e propô-la em um artigo científico.