pi_from_flips.R

# Adapted from Matt Parker's video here: https://www.youtube.com/watch?v=kahGSss6SsU
# Based on this paper here: https://arxiv.org/abs/2602.14487

library(tidyverse)
library(readxl)

# A "run" ends when cumulative heads first exceed 50% of flips
# The expected value of heads/total at that stopping time is pi/4

# Extract the heads/total ratio at the end of each run
extract_ratios <- function(flips) {
  n <- length(flips)
  ratios <- numeric(n)
  n_ratios <- 0L
  heads <- 0L
  total <- 0L

  for (flip in flips) {
    heads <- heads + flip
    total <- total + 1L
    if (heads * 2L > total) {
      n_ratios <- n_ratios + 1L
      ratios[n_ratios] <- heads / total
      heads <- 0L
      total <- 0L
    }
  }

  ratios[seq_len(n_ratios)]
}

estimate_pi <- function(flips) mean(extract_ratios(flips)) * 4

# Here are Matt Parker's original 10,000 flips
flips_excel <- read_excel(
  "~/Downloads/coinflip-stats.xlsx",
  sheet = "10k FLIPS",
  col_names = FALSE,
  col_types = "text",
  skip = 3
) |>
  select(1:3) |> 
  setNames(c("time", "started", "landed")) |> 
  drop_na() |> 
  filter(landed != "e") |> 
  mutate(landed = as.integer(landed == "h")) |> 
  pull(landed)

# And here's almost-pi using just these flips!
estimate_pi(flips_excel)
#> [1] 3.226593


# Here's pi from 10 million random coin flips
withr::with_seed(123456, {
  rbinom(10000000, 1, 0.5) |> estimate_pi()
})
#> [1] 3.131381