St. Powerball Paradox

I’ve been somewhat swept away with “lottery fever” the past few weeks. It hits at my house doubly, because we’re all interested in dreaming about “hitting it big,” but it’s professionally interesting for me as well. I always seem to find myself at the intersection of human behavior and advanced mathematics.

I bought a lottery ticket. Many of my co-workers think that I am opposed to playing the lottery, but this is not true. The pleasure that we derived from our one ticket was “on par” with other things (I replaced my wife’s weekly flowers with lotto tickets at her request). It also led to a fantastic discussion with my daughter and chief assistant of three “extensions” to the lottery. Fast forward:

Lottery tickets -> sequential lottery -> St. Petersburg Paradox -> How many heads in a row can RNG (random number generation) in R produce?

I was in the middle of discussing the St. Petersburg Paradox with my daughter and thought it might be interesting to demonstrate it with computer simulation, but I awoke in the middle of the night realizing that I couldn’t!

Warning: The rest of this article is about something we use a lot but don’t think about very often – RNGs are the “gateway drug” to the forbidden pleasures of Number Theory. It can lead to late nights huddled over old texts and contemplating 2ⁿ-1.

What I initially thought of doing was this: I would code up the St. Petersburg Paradox, run it a few million times, and then see how I did “on average.” Flipping a coin is random enough for me to call it random: small, unpredictable changes in air temperature, pressure and currents. There is no physical limit to the number of consecutive heads, and we expect the next flip to be “heads” with a probability of 50 percent regardless of the history of heads to date.

Using a RNG is not like flipping a coin, and the RNG is not random! The default RNG in R, which is a pretty good one, is the Mersenne Twist and is statistically random by virtue of passing randomness tests, known as the “Diehards.”

We can convert the uniform random number that lives on (0,1) by interpreting a “head” as runif(1) > .5.

We know that there is a limit to the longest sequence that can be in R; several estimates are possible, staying within our five-minute allocation, we know for certain that it is less than (see following equation). The full cycle of the RNG in R is unimaginably huge. Consider this:

The estimated time until the universe suffers “heat death” under closed universe theories is 10⁹⁰ years from now. Iterating the RNG at one per nanosecond until the end of time won’t get you the full cycle – not even close. (I was tempted to see if there was enough silicon in the universe to build enough computers to do the processing in parallel, but it is too depressing to think about.) One interesting feature of R is that the first time the RNG is invoked it pulls the seed from the CPU Clock, so almost all of us in the world are running a different instance of the R RNG with non-overlapping streams. In my experiments, looking for runs of “heads” a billion (a paltry ?2³⁰) or so at a time are laughably small (Figure 1).

The longest stream of “heads” I was able to produce in R was 33, which equates to winnings of “only” $537 million – which is large, but not infinite.

I will give bragging rights, and a bottle of Laphroig Cask Strength (provided it is legal and ethical to do so) for the reply that most substantially improves on my 33-ple that I found in R without using degenerate seeds (see below). They will need to send me both the string of random numbers and the .Random.seed vector for verification.

A few things to note:

1. Very important: When you start to manipulate .Random.seed, you are meddling with powers you cannot possibly comprehend. I strongly recommend beginning with safeseed = .Random.seed to store the current value, so that your RNG doesn’t become degenerate. Should you mess up your RNG, you can use .Random.seed = safeseed to restore operations. (Note: this is a good rule of practice when playing with the inner settings of any computer.)

2. Speaking of degeneracy, a vector of NAs will produce a long string of .5314 in R (see below). Using the command .Random.seed = as.integer(c(403, 444, rep(.Machine$integer.max-1, 624) will produce a long sequence of .9333. I’m defining these as degenerate because I think (without proof) that they are unreachable from a “normal” starting state. This doesn’t win a bottle.

The second point here is interesting. When I first started playing with this, I thought that the string of constant output was unending. It is not. The RNG eventually recovers (remember that the seed length is a 626-tuple, not a single “long” like many others. Consider the three graphs shown in Figure 2. 

A “real” number theorist would derisively call this an exercise in “number factualism.” I don’t necessarily disagree. I am not aware of any theory that could be brought against the problem of “how long will an RNG stay in a particular region,” but of course, this doesn’t mean it doesn’t exist.