Introduction


Nobody can be lucky all the time;
don’t think you’ve been abandoned in your prime,
but rather that you’re saving up your ration.

Piet Hein, Grooks


This markdown will simulate coin tosses and ask simple questions about runs of consectuive heads or tails. We will start by asking if luck predicts more luck, which will allow us to discover the concept of the Gambler’s Fallacy. We will then look at maximum runs of luck given a specific number of coin tosses, which will allow us to come across the concept of the fat-tailed distribution, which in turn is very common in nature.

Runs of luck 1: does luck predict more luck?

This will be the first section where we take the simulated coin tosses and analyze instances where there are a run of heads or a run of tails. We’re going to attack the concept of the Gambler’s Fallacy first. The problem and its magnitude is explained below:


It’s called the ‘hot hand’ fallacy – a belief that your luck comes in streaks – and it can lose you a lot of money. Win on roulette and your chances of winning again aren’t more or less – they stay exactly the same. But something in human psychology resists this fact, and people often place money on the premise that streaks of luck will continue – the so called ‘hot hand’.

The opposite superstition is to bet that a streak has to end, in the false belief that independent events of chance must somehow even out. This is known as the gambler’s fallacy, and achieved notoriety at the Casino de Monte-Carlo on 18 August 1913. The ball fell on black 26 times in a row, and as the streak lengthened gamblers lost millions betting on red, believing that the chances changed with the length of the run of blacks.

Why We Gamble Like Monkeys, BBC


You can see that Piet Hain’s poem this post opened up with falls into this category. Below, I wrote a function that flips a coin until there is a run of a specified number of heads in a row. The function flips a coin one more time after that and records only that last result. If 5 heads in a row leads to the 6th flip being also heads more often than not, then we will see it. If this is not the case, then there will be nearly as many heads as tails.

library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ ggplot2 3.3.5     ✓ purrr   0.3.4
## ✓ tibble  3.1.6     ✓ dplyr   1.0.8
## ✓ tidyr   1.2.0     ✓ stringr 1.4.0
## ✓ readr   2.1.2     ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
#' @title n run trial
#' @description Flips a fair coin (0 or 1 random number generator) over and over until the number of consective 1's reaches n times. Then flips a coin one more time. Returns a vector that contains the initial flips plus the final flip. 
#' @param n The number of consecutive 1's that will halt the iteration
#' @return A list where the first element is the initial binary vector generated and the second element is the flip that happens at the end of the specified run of n-heads consecutively generated.
NRunTrial <- function(n) {
    curr_run <- 0
    count <- 0
    curr_flips <- c()
    while (curr_run < n) {
        count <- count + 1
        flip <- sample(x = c(0, 1),
                   size = 1,
                   replace = FALSE)
        curr_flips <- c(curr_flips, flip)
        
        if(count > 1) {
            if (flip == 1 & curr_flips[count - 1] == 1) {
                curr_run <- curr_run + 1
            } else if(flip == 1 & curr_flips[count - 1] == 0) {
                curr_run <- 1
            } else if(flip == 0) {
                curr_run <- 0
            }
        } else {
            if(flip == 1) {
                curr_run <- 1
            } else {
                curr_run <- 0
            }
        }
    }
    flip <- sample(x = c(0, 1),
                   size = 1,
                   replace = FALSE)
    return(list(curr_flips, flip))
}

test <- NRunTrial(n = 5)
test[[1]]
##  [1] 0 0 0 1 1 0 0 1 0 0 1 1 1 1 0 0 0 1 0 0 1 1 1 1 0 0 1 1 1 0 1 0 0 1 0 0 1 1
## [39] 1 0 0 1 0 0 1 0 0 0 1 1 1 1 0 0 0 0 1 1 1 1 1
test[[2]]
## [1] 1

Now we run the function over and over. We use the same loop we made to produce our original bimodal distribution for the 100-toss Bernouli process. However, each toss in this case is actually the toss that occurrs after we’ve gotten 3 heads in a row. If this run of three heads “supercharges” the coin with luck, then the resulting distribution should be shifted or tailed to the right somehow, indicating that the 3 heads in a row leads to the next toss being heads more often.

set.seed(1)

n_sessions <- 100000

# Wait until we get three consecutive heads, then flip again. 
num_heads <- lapply(seq(n_sessions), function(i) {
    result <- lapply(seq(100), function(j) {
        NRunTrial(n = 3)[[2]]
    }) %>% unlist()
    result <- sum(result)/length(result)
    return(result)
}) %>% unlist()

# Histogram
hist(num_heads, breaks = 50, xlim = c(0, 1))

mean(num_heads)
## [1] 0.4996731
sd(num_heads)
## [1] 0.04993156

We see that we get the same bimodal distribution as we got before with a mean of 50% and a standard deivation of 5%. This means that there is no luck from previous runs being transferred to the next run. This result is an example of a bigger concept called the Gambler’s Fallacy. It’s the idea that a run of luck in a casino game doesn’t increase the odds that you’ll get lucky again.

If gamblers ran these computational experiments prior to going to the casino, perhaps they would save a little bit of money.

Runs of luck 2: maximum runs


“I was one of the luckiest NBA players ever to play with Michael Jordan, Tim Duncan, David Robinson and some of the greatest players ever,” Warriors coach Steve Kerr said. “As many spectacular things as Michael did, which he did nightly, I never saw him do that.” Steve Kerr on NBA star Klay Thompson’s record-setting 37 point quarter that stunned the basketball world on Jan. 23, 2015. —

When flipping coins, sometimes you get “runs of luck” where you have a string of heads or a string of tails that seems to defy the odds. We see this often in basketball when a certain player is “on fire” and can’t seem to miss a basket. Back to the coin tosses, we can calculate the probability of getting n heads in a row as multiplying 1/2 n times, or:

\[ 1/2^{n} \]

We’re going to ask a slightly different but related question here: if we flip a coin n times, what is the maximum run of luck we will get? Note that framing the question this way allows to examine the distribution of runs of luck given a specific number of tosses. This should give you some intuition around what the odds are that you’ll defy the odds.

We program the run of luck counting below, to initially get at the relationship between number of flips and maximum run of luck, which should be similar to 1/2^n probability calculation we introduced above.

set.seed(1)

#' @title Find Maximum Luck
#' @description Given a vector of coin flips (binary choices), find the largest number of 1's of 0's in a row
#' @param flips a numeric binary vector representing coin flips
#' @return a single integer denoting the maximum number of 1 or 0 in a row
FindMaxLuck <- function(flips) {
    previous_flip <- 0
    curr_count <- 0
    max_count <- 0
    
    for(i in flips) {
        if(i == previous_flip) {
            curr_count <- curr_count + 1
            if(curr_count > max_count) {
                max_count <- curr_count
            }
        } else {
            previous_flip <- i
            curr_count <- 1
        }
    }
    return(max_count)
}

test <- c(0, 1, 1, 0, 0, 0, 1, 0, 0)
FindMaxLuck(test)
## [1] 3
flip_number <- c(100, 1000, 10000, 100000, 1000000) 
flip_tib <- lapply(flip_number, function(i) {
    curr_flips <- sample(x = c(0, 1), size = i, replace = TRUE)
    result <- list(num_flips = i, max_luck = FindMaxLuck(curr_flips))
    return(result)
}) %>% bind_rows()


flip_tib
## # A tibble: 5 × 2
##   num_flips max_luck
##       <dbl>    <dbl>
## 1       100        7
## 2      1000       10
## 3     10000       16
## 4    100000       14
## 5   1000000       20

We observe that as we go up by orders of magnitude, the maximum luck only seems to increase by 3-5 at a time. If we take the number 20 at the bottom of the table, and ask what is the probability of getting 20 flips in a row, the answer is:

1/2^20
## [1] 9.536743e-07

9 in 10 million, which is close to 10 in 10 million, which is 1 in 1 million, which checks out on the table above. This is a good sanity check for our simulation experiment.

We’re now going to derive a “runs of luck” distribution, the purpose of framing our original question in the way we framed it. We’ll hypothesize that this distribution looks something like the binomial distribution, in that it’s a symmetric bell curve (spoiler alert: don’t hold your breath). Below, we simulate Bernouli processes of various sizes from 100 flips to several orders of magnitude above that. For each Bernouli process, we record the longest run of heads OR tails. We make no distinction between whether the longest run was heads or tails. We just care about the length of the longest run of either result.

flip_number <- c(100, 1000, 10000, 100000, 1000000) 
flip_list <- lapply(flip_number, function(i) {
    
    flips <- lapply(seq(10000), function(j) {
        result <- sample(x = c(0, 1), size = i, replace = TRUE) %>% FindMaxLuck()
        return(result)
    }) %>% unlist()
    
    return(flips)
    
}) 
names(flip_list) <- flip_number

count <- 0
invisible(lapply(flip_list, function(i) {
    count <<- count + 1
    hist(i, breaks = 20, 
         main = paste("max luck in", flip_number[count], "flips"), 
         xlim = c(0, 30))
}))

We can see that we were right in that it does form a bell curve, but we were wrong in that it is not symmetric: it has a longer tail to the right. We can see that as the number of flips increases, the distribution shifts to the right and seems to maintain it’s shape. We output summary statistics to get a better sense of how much it shift.

These distributions are known as fat-tailed distributions, which are relatively common in nature, and they have major implications in fields like risk analysis. In general, this fat-tailed distribution belongs to a family of distributions called extreme value distributions, which is part of a branch of math called extreme value theory. This is the type of math used to predict extreme events, like 100 year floods, or epidemics. The extreme value distribution might show up if you look at the longest run of made free throws for a given basketball player, but I’ll leave that as an exercise for the reader. All this being said, extreme value theory is very important branch of mathematics, which we ran head-first into by simply looking at runs of luck in our coin toss simulations.

We originally observed that for \(10^{n}\) tosses, the maximum run of luck seemed to be a one or two more than 3n. Below, we can generate the summary statistics, and then look at this relation with more precision.

flip_tib <- lapply(names(flip_list), function(i) {
    curr_flips <- flip_list[[i]]
    desc_stats <- summary(curr_flips) %>% as.list()
    result <- c(num_flips = i, sd = sd(curr_flips), desc_stats)
    return(result)
}) %>% bind_rows()

flip_tib$num_flips <- as.numeric(flip_tib$num_flips)
flip_tib
## # A tibble: 5 × 8
##   num_flips    sd  Min. `1st Qu.` Median  Mean `3rd Qu.`  Max.
##       <dbl> <dbl> <dbl>     <dbl>  <dbl> <dbl>     <dbl> <dbl>
## 1       100  1.79     3         6      7  6.99         8    21
## 2      1000  1.87     6         9     10 10.3         11    24
## 3     10000  1.86     9        12     13 13.6         15    26
## 4    100000  1.86    13        16     17 16.9         18    29
## 5   1000000  1.89    16        19     20 20.3         21    32
# Check relationship
qplot(log10(flip_tib$num_flips), flip_tib$Mean) + labs(x = "Number of zeros of flips", y = "Mean longest run of luck")

# Check for log10 linearity (run is always 1)
mean_slope <- lapply(seq(nrow(flip_tib)), function(i) {
    if(i == 1) {
        return(0)
    }
    result <- flip_tib[i,]$Mean - flip_tib[i - 1,]$Mean
    return(result)
}) %>% unlist()

mean_slope[-1]
## [1] 3.3202 3.3091 3.3112 3.3195
plot(mean_slope[-1])

The relationship between the number of zeros and the mean longest run of luck is linear. Note here the team “number of zeros” is a log10 transformation of a given number, as you see in the code. Number of zeros is simpler to write, so we are sticking with this terminology. Based on the second graph, the slope appears to be converging to between 3.3 and 3.4. In other words, it is not linear to begin with, but it appears to be converging to a line with a slope of around 3.35. For large numbers, we set our equation to y = 3.35x. This means if you flip a coin one million times (which we did), then the average maximum run of luck will be 6*3.35 = 20.1. If we check empirically, the answer is close, at 20.3. For a given project, you might need more precision than this, or you might simply need to know it’s the ballpark of 3 times the number of zeros plus a few, which we were able to quickly derive before we even were examining the extreme value distributions. Furthermore, the distributions we arrived at here can tell us much more information, like the minimum and maximum run of luck one would expect given a certain number of coin tosses.

To wrap up this section, we started by asking a simple question about maximum runs of luck given a number of coin tosses. Doing the calculations with simple code led us to a simple almost linear relation between the number of zeros in the coin tosses and the maximum run of luck, which led us to the often used log transformation of data, which then led us to the extreme value distribution and extreme value theory, which introduced us to the concept of fat-tailed distributions, which are very common in nature. This gives us some potential intuition around extreme value events, like runs of luck in shooting a basketball. All of this from asking a simple question about our coin toss simulations.