How to simulate the St. Petersburg paradox 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))


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.
