Simulation of Skulls and Roses bluffing game

skull.py

import numpy as np
from numba import jit, int32
from numpy.random import default_rng

rng = default_rng()


def startHands(numplayers, cards=4):
    return np.stack((cards * np.ones((numplayers,)), np.ones((numplayers,))))


# If a player has N cards, with 1 being a skull and he places P cards down
# The chance that 1 is a skull is:
# 1 - \prod_{i=0}^{P-1} (N-i-1)/(N-i)
# This simplifies to: P/N


def playCards(samples, hands, nCards):
    nPlayers = hands.shape[1]
    out = np.zeros(nCards)

    nCardsEach = nCards / nPlayers
    frac = nCardsEach - int(nCardsEach)
    cardsPerPlayer = int(nCardsEach) + (np.arange(nPlayers) / nPlayers < frac)

    skullsEach = rng.binomial(
        1, cardsPerPlayer / hands[0, :] * hands[1, :], (samples, nPlayers)
    )
    skullsTotal = np.sum(skullsEach, axis=1)
    rosesTotal = nCards - skullsTotal

    return skullsTotal, rosesTotal


def draw(samples, skulls, roses, bet):
    skullsDrawn = rng.hypergeometric(skulls, roses, bet, (skulls.shape[0], samples))
    return skullsDrawn

# Returns the probability of winning the round
def prob(hands, place, bet=None, N=1000):
    if bet is None: bet = place
    skulls, roses = playCards(N, hands, place)
    drawn = draw(N, skulls, roses, bet)
    return 1 - np.count_nonzero(drawn) / (N*N)

N = 1000
P = 4
C = 4
D = C # 4
prob(startHands(P), C, D, N)