Binomial distribution simulation python

Question

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()}

OK, look at my code. Your code has obvious bug in np.random.uniform(size=(nsims, 7)) to begin with... — Severin Pappadeux
– Severin Pappadeux, Commented Sep 19, 2020 at 16:12
You're welcome. If you're satisfied with the answer you could accept and/or bump it. — Severin Pappadeux
– Severin Pappadeux, Commented Sep 20, 2020 at 20:43
And I updated the code with simple example how to compute expected number of games exactly and via Monte Carlo — Severin Pappadeux
– Severin Pappadeux, Commented Sep 21, 2020 at 3:30

Severin Pappadeux · Accepted Answer · 2020-09-21 03:29:52Z

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:

Get exact answers
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))

Collectives™ on Stack Overflow

Binomial distribution simulation python

1 Answer 1

Comments

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

Comments

Related