plot.py

import numpy as np
import matplotlib.pyplot as plt
from math import comb

def tax_coefficient(r):
    if r <= 100:
        return 0
    elif 100 < r <= 200:
        return 0.025
    elif 200 < r <= 500:
        return 0.05
    elif 500 < r <= 2500:
        return 0.1
    else:
        return 0.2

def compute_ev(N, JACKPOT):
    p_51 = (1 / comb(45, 5)) * (1 / 20)
    r_51 = JACKPOT

    p_5 = (1 / comb(45, 5)) * (19 / 20)
    r_5 = 100_000

    p_41 = (comb(5, 4) * comb(40, 1) / comb(45, 5)) * (1 / 20)
    r_41 = 2_500

    p_4 = (comb(5, 4) * comb(40, 1) / comb(45, 5)) * (19 / 20)
    r_4 = 50

    p_31 = (comb(5, 3) * comb(40, 2) / comb(45, 5)) * (1 / 20)
    r_31 = 50

    p_3 = (comb(5, 3) * comb(40, 2) / comb(45, 5)) * (19 / 20)
    r_3 = 2

    p_21 = (comb(5, 2) * comb(40, 3) / comb(45, 5)) * (1 / 20)
    r_21 = 2

    p_11 = (comb(5, 1) * comb(40, 4) / comb(45, 5)) * (1 / 20)
    r_11 = 1.5

    p_2 = (comb(5, 2) * comb(40, 3) / comb(45, 5)) * (19 / 20)
    r_2 = 1

    expected_n_51_winners = N * p_51
    expected_n_5_winners = N * p_5

    ev = (
        p_51 * (r_51 / (1+expected_n_51_winners)) * (1 - tax_coefficient(r_51 / (1+expected_n_51_winners))) +
        p_5 * (min(r_5, 2_000_000 / (1+expected_n_5_winners))) * (1 - tax_coefficient(min(r_5, 2_000_000 / (1+expected_n_5_winners)))) +
        p_41 * r_41 * (1 - tax_coefficient(r_41)) +
        p_4 * r_4 * (1 - tax_coefficient(r_4)) +
        p_31 * r_31 * (1 - tax_coefficient(r_31)) + 
        p_3 * r_3 * (1 - tax_coefficient(r_3)) + 
        p_21 * r_21 * (1 - tax_coefficient(r_21)) + 
        p_11 * r_11 * (1 - tax_coefficient(r_11)) + 
        p_2 * r_2 * (1 - tax_coefficient(r_2))
    ) - 1
    return ev

# Define ranges for N and JACKPOT
N_vals = np.linspace(1_000_000, 20_000_000, 30)
JACKPOT_vals = np.linspace(1_000_000, 100_000_000, 30)
N_grid, JACKPOT_grid = np.meshgrid(N_vals, JACKPOT_vals)
ev_grid = np.vectorize(compute_ev)(N_grid, JACKPOT_grid)

# Plot
fig = plt.figure(figsize=(10,7))
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(N_grid, JACKPOT_grid, ev_grid, cmap='viridis')
ax.set_xlabel('N (Coupons Played)')
ax.set_ylabel('JACKPOT')
ax.set_zlabel('EV')
ax.set_title('EV as a function of N and JACKPOT')
plt.tight_layout()
plt.show()

tzoker.py

from math import comb

N = 7_000_000

def tax_coefficient(r):
    if r <= 100:
        return 0
    elif 100 < r <= 200:
        return 0.025
    elif 200 < r <= 500:
        return 0.05
    elif 500 < r <= 2500:
        return 0.1
    else:
        return 0.2

# p_{event} = probability of event
# r_{event} = reward of event

p_51 = (1 / comb(45, 5)) * (1 / 20)
r_51 = 22_000_000

p_5 = (1 / comb(45, 5)) * (19 / 20)
r_5 = 100_000

p_41 = (comb(5, 4) * comb(40, 1) / comb(45, 5)) * (1 / 20)
r_41 = 2_500

p_4 = (comb(5, 4) * comb(40, 1) / comb(45, 5)) * (19 / 20)
r_4 = 50

p_31 = (comb(5, 3) * comb(40, 2) / comb(45, 5)) * (1 / 20)
r_31 = 50

p_3 = (comb(5, 3) * comb(40, 2) / comb(45, 5)) * (19 / 20)
r_3 = 2

p_21 = (comb(5, 2) * comb(40, 3) / comb(45, 5)) * (1 / 20)
r_21 = 2

p_11 = (comb(5, 1) * comb(40, 4) / comb(45, 5)) * (1 / 20)
r_11 = 1.5

p_2 = (comb(5, 2) * comb(40, 3) / comb(45, 5)) * (19 / 20)
r_2 = 1

print(f"p_51 = {p_51}")
print(f"p_5 = {p_5}")
print(f"p_41 = {p_41}")
print(f"p_4 = {p_4}")
print(f"p_3 = {p_3}")
print(f"p_51 = {p_51}")
print(f"p_51 = {p_51}")
print(f"p_51 = {p_51}")
print(f"p_51 = {p_51}")

# Only relevant for top 2 categories
# expected # winners  of a prize category = E[Binomial(n, p)] = n * p
expected_n_51_winners = N * p_51
expected_n_5_winners = N * p_5

ev = (
    p_51 * (r_51 / (1+expected_n_51_winners)) * (1 - tax_coefficient(r_51 / (1+expected_n_51_winners))) +
    p_5 * (min(r_5, 2_000_000 / (1+expected_n_5_winners))) * (1 - tax_coefficient(min(r_5, 2_000_000 / (1+expected_n_5_winners)))) +
    p_41 * r_41 * (1 - tax_coefficient(r_41)) +
    p_4 * r_4 * (1 - tax_coefficient(r_4)) +
    p_31 * r_31 * (1 - tax_coefficient(r_31)) + 
    p_3 * r_3 * (1 - tax_coefficient(r_3)) + 
    p_21 * r_21 * (1 - tax_coefficient(r_21)) + 
    p_11 * r_11 * (1 - tax_coefficient(r_11)) + 
    p_2 * r_2 * (1 - tax_coefficient(r_2))
) - 1

print("ev = ", ev)