Search code examples
pythonrandomprobabilitymontecarlo

Binomial distribution simulation python


Suppose 2 teams A and B are playing a series of games and the first team to win 4 games wins the series. Suppose that team A has a 55% chance of winning each game and that the outcome of each game is independent.

(a) What is the probability that team A wins the series? Give an exact result and confirm it via simulation.

(b) What is the expected number of games played? Give an exact result and confirm it via simulation.

(c) What is the expected number of games played given that team A wins the series? Give an exact result and confirm it via simulation.

(d) Now suppose we only know that team A is more likely to win each game, but do not know the exact probability. If the most likely number of games played is 5, what does this imply about the probability that team A wins each game?

This is what I have done but not getting it..need some input. Thank you

import numpy as np

probs = np.array([.55 ,.45])
nsims = 500000

chance = np.random.uniform(size=(nsims, 7))

teamAWins = (chance > probs[None, :]).astype('i4')
teamBWins = 1 - teamAWins

teamAwincount = {}
teamBwincount = {}
for ngames in range(4, 8):
    afilt = teamAWins[:, :ngames].sum(axis=1) == 4
    bfilt = teamBWins[:, :ngames].sum(axis=1) == 4

    teamAwincount[ngames] = afilt.sum()
    teamBwincount[ngames] = bfilt.sum()

    teamAWins = teamAWins[~afilt]
    teamBWins = teamBWins[~bfilt]

teamAwinprops = {k : 1. * count/nsims for k, count in teamAwincount.iteritems()}
teamBwinprops = {k : 1. * count/nsims for k, count in teamBwincount.iteritems()}

Solution

  • Ok, here is the idea and the code to get you going.

    This is I believe is Negative Binomial Distribution, which is quite easy to implement, and compute probabilities for favorite and underdog.

    With that code you have whole set of events defined, with probabilities which properly sum to 1. From that you could:

    1. Get exact answers
    2. Check simulation against probabilities

    Simulation code added the counters for many events, and single event simulator. So far it looks like it is making about the same probabilities as Negative binomial formula

    Code, Python 3.8 x64 Win10

    import numpy as np
    import scipy.special
    
    # Negative Binomial as defined in
    # https://mathworld.wolfram.com/NegativeBinomialDistribution.html
    def P(x, r, p):
        return scipy.special.comb(x+r-1, r-1)*p**r*(1.0-p)**x
    
    def single_event(p, rng):
        """
        Simulates single up-to-4-wins event,
        returns who won and how many opponent got
        """
        f = 0
        u = 0
        while True:
            if rng.random() < p:
                f += 1
                if f == 4:
                    return (True,u) # favorite won
            else:
                u += 1
                if u == 4:
                    return (False,f) # underdog won
    
    def sample(p, rng, N):
        """
        Simulate N events and count all possible outcomes
        """
    
        f = np.array([0, 0, 0, 0], dtype=np.float64) # favorite counter
        u = np.array([0, 0, 0, 0], dtype=np.float64) # underdog counter
    
        for _ in range(0, N):
            w, i = single_event(p, rng)
            if w:
                f[i] += 1
            else:
                u[i] += 1
    
        return (f/float(N), u/float(N)) # normalization
    
    def expected_nof_games(p, rng, N):
        """
        Simulate N events and computes expected number of games
        """
    
        Ng = 0
        for _ in range(0, N):
            w, i = single_event(p, rng)
    
            Ng += 4 + i # 4 games won by winner and i by loser
    
        return float(Ng)/float(N)
    
    
    p = 0.55
    
    # favorite
    p04 = P(0, 4, p)
    p14 = P(1, 4, p)
    p24 = P(2, 4, p)
    p34 = P(3, 4, p)
    
    print(p04, p14, p24, p34, p04+p14+p24+p34)
    
    # underdog
    x04 = P(0, 4, 1.0-p)
    x14 = P(1, 4, 1.0-p)
    x24 = P(2, 4, 1.0-p)
    x34 = P(3, 4, 1.0-p)
    print(x04, x14, x24, x34, x04+x14+x24+x34)
    
    # total probability
    print(p04+p14+p24+p34+x04+x14+x24+x34)
    
    # simulation of the games
    rng = np.random.default_rng()
    f, u = sample(p, rng, 200000)
    print(f)
    print(u)
    
    # compute expected number of games
    
    print("expected number of games")
    E_ng = 4*p04 + 5*p14 + 6*p24 + 7*p34 + 4*x04 + 5*x14 + 6*x24 + 7*x34
    print(E_ng)
    # same result from Monte Carlo
    print(expected_nof_games(p, rng, 200000))