Binomial distribution simulation python

user_14217438 0
Viewed 968

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()}
montecarlo probability random python
1 Answers

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))
user_40446960
Related