Previously in this series:
 The “lost boarding pass” puzzle
 The “deadly board game” puzzle
 The “knight on an infinite chessboard” puzzle
 The “largest stock profit or loss” puzzle
 The “birthday paradox” puzzle
I have an interest in probability puzzles and riddles, and especially in simulating them in R. I recently learned about Feller’s cointossing puzzle, from the book Mathematical Constants by Steven Finch. (I recommend the book if you like the topic too!)
Mathematician William Feller posed the following problem:
If you flip a coin times, what is the probability there are no streaks of heads in a row?
Note that while the number of heads in a sequence is governed by the binomial distribution, the presence of consecutive heads is a bit more complicated, because the presence of a streak at various points in the sequence isn’t independent. This reminds me a bit of one of my earlier tidyverse simulations:
A #tidyverse simulation to demonstrate that if you wait for two heads in a row, it takes 6 flips on average, while you wait for a heads then a tails, it takes 4 flips on average
— David Robinson (@drob) June 17, 2018
h/t @CutTheKnotMath #rstats pic.twitter.com/V0zgOmCy7t
To continue my series of simulating probability puzzles in the tidyverse, I’d like to show how we’d approach simulating Feller’s cointossing problem, and comparing it to the exact values. (In the process, we also see how we’d calculate a Fibonacci sequence in one line!)
Simulating a single sequence
Let’s start with values : what’s the probability that a sequence of 20 flips contains no streaks of length 3? You can flip a sequence of coins with rbinom()
.
In this case, there were indeed a few streaks of 3 heads in a row. How could determine that in R?
Well, we could use dplyr’s window function lead()
(which moves each flip forward one in the sequence), to ask if there are any flips sets in which a coin, the next coin, and the one after that are all 1 (heads).
Indeed, there are (though notice the last two are NA, since there is no lead()
coin).
Remember that Feller was looking for the probability there are no streaks in the sequence. We use !any()
to check this:
This gives us an approach that, similar to our previous tidyverse approaches to simulation, we can repeat and summarize across parameter values using tidyr’s crossing()
and purrr’s map_lgl()
.
This takes about 5 seconds on my machine. We can then visualize the probability that there are no streaks as a function of the number of flips in the sequence.
It looks like for , the probability there are no sequences of three is 7/8 for a sequence of 3 flips, crosses 50% roughly when there are 10 flips, and then is rather close to zero by the time there are 50 flips.
Extending for multiple values of k
Once we’re not fixed to , we can’t use x & lead(x) & lead(x, 2)
to check for the presence of a streak anymore.^{1} As a replacement, I’d like to introduce a useful base R function called rle
, for “runlength encoding”.
A runlength encoding divides a vector down into streaks of consecutive values. It turns the vector into two components: the lengths of each streak, and the value in each. We can use these in combination !any(r$values & r$lengths >= len)
 to check if there are any streaks of heads greater than a certain length. (This is a good example of how knowing slightly obscure base R functions, like rle
, gives you a toolbox for elegant and efficient solutions).
By adding a value k
to our crossing()
, we can then visualize the probability for each value of k.
The longer the streak, the less likely the sequence won’t contain it, which makes sense. By the time the sequence is length 40, it’s almost certain to contain a stretch of 2 heads, very likely to contain a stretch of 3 heads, and has a 75% chance to contain a stretch of 4 heads.
Feller’s cointossing constants
Something I like about simulations is that they can doublecheck mathematical results.
When Feller looked at this problem, he proved a result about , namely:
Where and are Feller’s constants. (You can find a few such values, and some more details, on Wikipedia). We could compare those exact values to the simulation, by creating a table of the constants and joining them.
Calculating the probability of Fibonacci numbers
But an approach I like even more than Feller’s constants is to calculate the exact probability based on higherorder Fibonacci sequences.
where is the term of the th order Fibonacci sequence. The first few paragraphs of this paper points out why. (In short for : the number of sequences of length that have no streaks of 2 is all the sequences of length that are followed by a , plus all the sequences of length that are followed by a . This is divided by the possible sequences.)
Let’s talk about Fibonacci sequences! Each step in a Fibonacci sequence is the sum of the previous 2, after starting with (1, 1). To get that in R, you’d keep applying the step c(., sum(tail(., 2)))
again and again (tail()
gets the last items of a vector).
This can be done in one line (trick of the day!) with the reduce
function from purrr, which calls a function for each element in a vector while passing along the result. When passed a dummy vector, like 1:50
, and an initial value, like c(1, 1)
(the first two), it’s a quick way to say “call this function 50 times”.
In higher order Fibonacci sequences, the terms are the sum of 3 (“tribonacci”), 4 (“tetranacci”), or more previous terms, meaning they grow even faster. We could create a function that calculates those series.
So returning to our simulation, we can confirm our math.
Notice what a wide range of tools can be used in simulations. Besides our usual collection of tidyverse tricks like crossing()
, we used rle()
(a handy trick any time you need to examine consecutive streaks), and reduce()
(useful for setting up recursive relationships like in the Fibonacci sequence).
I’m really enjoying these probability puzzle simulations. If you have a favorite probability puzzle you’d like me to simulate, please put in the comments!

With
reduce
, there actually is a way we could take thelead()
approach with an arbitrary streak length (left as an exercise to the reader!). But I found it’s about 10X slower than therle()
approach above, so I’m focusing on this one. ↩