Modeling The St. Petersburg Paradox

Background:

The St. Petersburg game is a game of chance; the person placing the wager is asked for a certain amount. After paying the casino to play a fair coin is flipped until a tails is reached. The player wins 2^{number of flips} dollars. The chances of winning the nth round are given by 1/2^n for example the chances of winning the eighth round are 1/256 and the payoff for that round is 256. The expected payoff for any particular round is the chance of winning multiplied by the payoff, for this game it is always 1. The expected payoff for a game is the sum of all possible rounds. In this case there are an infinite number of rounds so the expected payoff is infinity.  Hypothetically then, the perfectly rational person will pay any amount in order to play, but how realistic is this?

Setup:

I decided to write a simulation in Python (of course) to actually check what the average player should pay to play.

#!/usr/bin/python
'''stpetersburg.py -- A saint petersburg simulator.

Copyright 2011 Joseph Lewis <joehms22@gmail.com>

MIT License http://www.opensource.org/licenses/mit-license.php

'''
import random

def play():
    total_won = 0
    payoff_matrix = {}

    print("Iter:\tFlips:\tWinnings:\tWon?")

    for i in range(iterations):
        winnings = 0
        flips = flip(max_game)  #Flip the coin until tails.

        #Determine winnings based upon the wager multiplier and flips.
        if flips:
            winnings = multiplier**(flips)

        #Stats
        total_won += winnings
        if flips in payoff_matrix.keys():
            payoff_matrix[flips] = payoff_matrix[flips] + 1
        else:
            payoff_matrix[flips] = 1

        print("%i\t%i\t%i\t\t%s" % (i, flips, winnings, winnings >= initial_wager))

    #Do some basic stats.
    avg_won = float(total_won) / iterations
    print ("Avg. Winnings:\t%d" % (avg_won))
    print ("Flips x Frequency")
    for k in payoff_matrix.keys():
        print("%s x %s" % (k, payoff_matrix[k]))
    print ("Money Won:   %s" % (total_won))
    print ("Money Waged: %s" % (iterations * initial_wager))


def flip(max_flip):
    '''Flips a coin until max_flip is reached or the coin turns up
    tails.  Returns the number of flips that came up heads.

    If max_flip is -1 flipping will continue until a tails is reached.
    '''
    numflips = 0

    while numflips != max_flip:
        numflips += 1
        if random.randint(0,1): #0 = heads 1 = Tails:
            return numflips


    return numflips


def ask_int(question, default):
    '''Asks the user for an int printing out the default and question.
    '''
    num = raw_input(question + " (Default: %s)\n" % (str(default)))
    if num == '':
        return default
    return int(num)


if __name__ == "__main__":
    print __doc__ #Notify the user of the license (top of file)

    #Set up simulation
    iterations = ask_int("How many iterations?", 100000)
    initial_wager = ask_int("What is the initial wager?", 2)
    multiplier = ask_int("What is the multiplier per win?", 2)
    max_game = ask_int("What is the bound on the number of games -1 for infinity?", -1)

    #Start simulation
    play()

Running:
This is an output of the default settings (which can be tweaked to easily accommodate other versions of the problem).

Avg. Winnings:  18
Flips x Frequency
1 x 50078
2 x 25019
3 x 12474
4 x 6286
5 x 3085
6 x 1455
7 x 791
8 x 396
9 x 197
10 x 110
11 x 51
12 x 32
13 x 17
14 x 5
15 x 2
16 x 1
18 x 1
Money Won:   1858488
Money Waged: 200000

The output was fairly consistent the several times I have run it, the amount won is about ten times the money wagered, so therefore you should want to bet no more than about twenty dollars.  There you have it, a rational cutoff for playing the paradox :)