Personal code snippets of @tmasjc

Image by Mads Schmidt Rasmussen from unsplash.com

Minimal Bootstrap Theme by Zachary Betz

# Puzzle: Biased Coin Toss

Apr 21, 2019 #monte-carlo

From this paper “What’s Past Is Not Prologue” comes a quiz,

You are presented with two coins: one is fair, and the other has a 60% chance of coming up heads. Unfortunately, you don’t know which is which. How many flips would you need to perform in parallel on the two coins to give yourself a 95% chance of correctly identifying the biased one?

The following demonstrates 2 simulation methods to obtain the result.

library(tidyverse)
library(ggthemes)
library(discreteRV)
set.seed(1212)

# ggplot theme
old <- theme_set(theme_tufte() + theme(text = element_text(family = 'Menlo')))
# limit range from 50 to 200 flips
trials <- 50:200

# Method One --------------------------------------------------------------

# actually simulate the action of tossing 2 coins
method_one <- function(x) {
regular = rbinom(x, 1, .5)
biased = rbinom(x, 1, .6)
# return
sum(biased) - sum(regular)
}

# repeat action and bind into a data frame
res.one <- map_dbl(trials, ~ {
sum(replicate(3000, method_one(.x)) > 0) / 3000
}) %>% bind_cols(ind = trials, prob = .)

# first occurance
(min.flip.one <- filter(res.one, prob >= 0.95)[1, ])
## # A tibble: 1 x 2
##     ind  prob
##   <int> <dbl>
## 1   136 0.952
# vis result
res.one %>%
ggplot(aes(ind, prob)) +
geom_line() +
geom_vline(xintercept = min.flip.one$ind, lty = 4, col = "red") + geom_hline(yintercept = min.flip.one$prob, lty = 4, col = "red") +
annotate("text", x = min.flip.one$ind, y = .5, label = as.character(min.flip.one$ind)) +
scale_y_continuous(breaks = seq(.5, 1, .1), labels = scales::percent) # Method Two --------------------------------------------------------------

# we will reuse this function repeatedly to test different diff trials
method_two <- function(x) {
# how many tosses?
regular <- RV(1:x, dbinom(1:x, size = x, prob = .5))
biased <- RV(1:x, dbinom(1:x, size = x, prob = .6))

# joint pmf
j2 <- joint(regular, biased)

# probability here can be understood as proportion of combination of possible outcomes
tibble(prob = probs(j2), comb = outcomes(j2)) %>%
separate(comb, c("fair", "biased"), ",", convert = TRUE) %>%
# if biased lands more than unbiased coin, choose correctly
# if both lands the same, choose randomly
mutate(correct = case_when(
biased > fair ~ 1,
biased == fair ~ 0.5,
biased < fair ~ 0
)) %>%
# the likelihood that guessing correctly happens
mutate(likelihood = prob * correct) %>%
filter(correct > 0) %>%
summarise(conf = sum(likelihood)) %>%
pull(conf)
}

# loop thru trials
res.two <- map_dbl(trials, ~ method_two(.x)) %>% bind_cols(ind = trials, prob = .)

# minimum flip to obtain 95% confidence
(min.flip.two <- filter(res.two, prob >= 0.95)[1, ])
## # A tibble: 1 x 2
##     ind  prob
##   <int> <dbl>
## 1   134 0.950
# vis result
res.two %>%
ggplot(aes(ind, prob)) +
geom_line() +
geom_vline(xintercept = min.flip.two$ind, lty = 4, col = "red") + geom_hline(yintercept = min.flip.two$prob, lty = 4, col = "red") +
annotate("text", x = min.flip.two$ind, y = .5, label = as.character(min.flip.two$ind)) +
scale_y_continuous(breaks = seq(.5, 1, .1), labels = scales::percent) 