How to simulate the St. Petersburg paradox

Question

I am trying to find the expected value of the game in the St. Petersburg paradox. The game has infinite expected value, but in my simulation the payouts are about 15 on average, which is way too small.

Increasing iterations does increase the simulated mean but not much. Can anyone tell me or point to resources how to properly get such simulations to converge to their expected value, where events with very low probability are important? This is the code I used:

payout_list <- list()
B <- 1000000

for(i in 1:B){
  payout <- 1
  tails <- TRUE
  
  while(tails){
    payout <- 2*payout
    tails <- sample(c(T,F), size = 1)
  }
  payout_list <- append(payout_list, payout)
}

avg <- mean(as.numeric(payout_list))

Answer 1

This is more a story of instability than of infinity.

You can spot the inexistence of a mean of a random variable $X$ by examining how the empirical mean changes over time. When the mean exists, eventually the mean of a long series of independent draws will settle down to the mean of $X:$ that is what various Laws of Large Numbers assert.

Now, although no finite simulation can definitively establish whether a mean exists or not, simulations can be suggestive.

As Wikipedia describes it,

A casino offers a game of chance for a single player in which a fair coin is tossed at each stage. The initial stake begins at 2 dollars and is doubled every time heads appears. The first time tails appears, the game ends and the player wins whatever is in the pot.

This game is certain to terminate, because the chance that it lasts longer than $n$ tosses is only $2^{-n},$ which can be made smaller than any positive number, showing that the chance the game does not terminate is less than all positive numbers: it must be zero.

The number of tosses follows a Geometric distribution, allowing for efficient simulation of this game, as in this R implementation of a million independent iterations.

X <- 1 + rgeom(1e6, 1/2)

A histogram of these results -- which are the binary logarithms of the winnings -- gives us a picture of their relative frequences. The steady decrease from $1/2$ for surviving just one roll, to $1/4,$ to $1/8,$ and so on, is apparent, indicating this simulation is doing the right thing.

The second panel plots mean winnings after each iteration. They don't appear to be settling down. The infinite expectation is manifest in the occasional large jumps. Think about what a jump must mean. The leap from near 20 to almost 25 around iteration 400,000, for instance, indicates the total winnings must suddenly have increased by about $(25 - 20)\times 4\times 10^5 \approx 2^{11}.$ It took a sequence of 11 tails in a row to do that. Every once in a while, going on as long as the game could ever be played, there will be such jumps.

Others have used simulations in the same way to analyze other situations with undefined or infinite means, such as a Cauchy random variable.

Answer 2

I'm not sure you can numerically simulate a process involving infinity. Every individual instance of the game terminates with a finite value, so any finite number of plays will result in the mean value also being finite. No matter how many times you simulate the game in practice, your simulation mean will never reach the expected value of the game, which is infinite.

The underlying issue is that a numerical simulation cannot converge to infinity - numerical models produce numerical results, but infinity isn't a number. At best, your simulation can return a very large number, but that does not itself imply that the process actually converges to an infinite expected value.