St. Petersberg Paradox

stpetersberg.py

# St. Petersberg paradox simulation
# https://en.wikipedia.org/wiki/St._Petersburg_paradox

import numpy as np

rng = np.random.default_rng(seed=42)


def geometric():
    """Play one game of 'st_petersberg paradox' by simulating coin throws.
    Returns the number of heads until tail. (i.e. geometric distribution with p=0.5)"""
    combo = 0
    while rng.random() <= 0.5:
        combo += 1
    return combo

def st_petersberg_paradox(it):
    """Play {it} games of 'st_petersberg paradox' by simulating coin throws.
    Returns the average reward. (0 -> 0, x -> 2**x for x > 0)"""
    rew = 0
    for _ in range(it):
        combo = geometric()
        rew += 2 ** combo if combo else 0  # 0, 2, 4, 8, 16, ...
    return rew/it

def st_petersberg_paradox_batch(it, batch_throw):
    """Play {it} games of 'st_petersberg paradox' by batch coin throws using numpy.
    Returns the average reward. (0 -> 0, x -> 2**x for x > 0)"""
    rew = 0
    combo = 0
    multiplier = 1.
    i = batch_throw
    cit = it
    while cit > 0:
        if i == batch_throw:
            throw = rng.random(batch_throw) > 0.5
            i = 0
        if throw[i]:
            rew += 2 ** combo if combo else 0
            cit -= 1
            combo = 0
            multiplier = 1.
        else:
            multiplier *= 2
            combo += 1
        i+=1
    return rew/it

def st_petersberg_paradox_geometric_batch(it, batch_throw):
    """Play {it} games of 'st_petersberg paradox' using batch geometric distribution.
    Returns the average reward. (0 -> 0, x -> 2**x for x > 0)"""
    rew = 0
    cit = it
    while cit > 0:
        current_throw = min(batch_throw, cit)
        throw = rng.geometric(0.5,current_throw)-1
        reward = np.where(throw==0, 0, 2.**throw)
        rew += np.sum(reward)
        cit -= current_throw
    return rew/it


if __name__ == '__main__':
    """
    1 4.0
    10 4.6
    100 2.34
    1000 6.21
    10000 10.0732
    100000 8.7018
    1000000 11.879142
    10000000 11.950668
    100000000 14.1001756
    1000000000 17.130889618
    """
    for i in range(10):
        #print(10**i, st_petersberg_paradox(10**i))
        #print(10**i, st_petersberg_paradox_batch(10**i, 1000000)) 
        print(10**i, st_petersberg_paradox_geometric_batch(10**i, 1000000))