Calculating the Probability of Becoming Rich in the St. Petersburg Paradox Game Before Going Broke

gamble.py

from scipy.sparse import lil_matrix
import scipy.sparse as sp
import scipy.sparse.linalg as linalg
import numpy
import math
import matplotlib.pyplot as plt
import matplotlib.widgets as wg 


def NextBadget(badget, fee, headCount):
    return badget - fee + ((2 ** (headCount -1)) if 1 <= headCount else 0)

# A of "x = Ax + b"
def CalcCoeffMatrix(fee, goal):
    coeffMatrix =lil_matrix((goal, goal))
    for badget in range(0, goal):
        if(badget < fee):
            continue
        headCount = 0
        nextBadget = NextBadget(badget, fee, headCount)
        while(nextBadget < goal):
            coeff = 0
            if(fee <= nextBadget):
                coeff =  0.5 ** (headCount + 1)
                coeffMatrix[badget, nextBadget] = coeff
            headCount = headCount + 1
            nextBadget = NextBadget(badget, fee, headCount)
    return coeffMatrix


# b of "x = Ax + b"
def CalcBoundArray(fee, goal):
    arr = numpy.zeros(goal)
    for badget in range(fee, goal):
        overGoalHeadCount = math.ceil(math.log2(goal + fee - badget) + 1)
        arr[badget] = 1 / (2 ** overGoalHeadCount)
    return arr

def Run(fee, goal):
    coeffMatrix = CalcCoeffMatrix(fee, goal)
    #print(coeffMatrix)


    # E-A of "(E-A)x = b"
    matrix = sp.identity(goal) - coeffMatrix
    #print(matrix)

    arr = CalcBoundArray(fee, goal)
    result = linalg.spsolve(matrix, arr)
    return result


Initfee = 1
InitGoal = 5
result = Run(Initfee, InitGoal)
print("result:", result)


fig = plt.figure(figsize=(12, 8))
ax = fig.add_axes([0.1,0.3,0.8,0.6])
ax.set_ylim(0, 1)
ax.set_xlim(0, len(result))
plot, = ax.plot(result)
ax_a = fig.add_axes([0.15, 0.1, 0.6, 0.04])
ax_p = fig.add_axes([0.15, 0.05, 0.6, 0.04])


Maxfee = 2 ** 11
sli_a = wg.Slider(ax_a, 'fee', 1, math.floor(math.log2(Maxfee)), valinit=Initfee,valstep=1)
sli_p = wg.Slider(ax_p, 'goal', 1, Maxfee, valinit=InitGoal,valstep=1)

def update(val):
   result = Run(sli_a.val, sli_p.val)
   plot.set_data(range(0, len(result)), result)
   ax.set_xlim(0, len(result))
   

sli_a.on_changed(update)
sli_p.on_changed(update)


plt.ylim(0,1)
plt.show()