monte carlo simulation python

Question

I would like to simulate a seven game baseball playoff series. Let's say I have the the win probabilities for each game in the series. I would like to know the probabilities for each possible series outcome. ie TeamA in 4 games, TeamB in 4 games, TeamA in 5 games, etc.

This is what I came up with and it seems to work but I think it could be done better.

winPercGM1 = .5
winPercGM2 = .56
winPercGM3 = .47
winPercGM4 = .55
winPercGM5 = .59
winPercGM6 = .59
winPercGM7 = .38
winPercs = [winPercGM1, winPercGM2,  winPercGM3, winPercGM4, winPercGM5, winPercGM6,     winPercGM7]

def WinSeries():   
    teamAwins = 0
    teamBwins = 0  
    for perc in winPercs:
        if teamAwins == 4:            
            break            
        elif teamBwins == 4:            
            break            
        elif perc > np.random.random():
            teamAwins += 1            
        else:
            teamBwins += 1            
    return teamAwins, teamBwins

def RunFun(n):
teamAWins = []
teamBWins = []   
for i in xrange(n):
    result = WinSeries()
    teamAWin = result[0]
    teamBWin = result[1]        
    teamAWins.append(teamAWin)
    teamBWins.append(teamBWin)       
return teamAWins, teamBWins


n = 500000
results = RunFun(n)

teamAwinSeries = results[0]
teamBwinSeries = results[1]

teamBin4 = teamAwinSeries.count(0)/n
teamBin5  = teamAwinSeries.count(1)/n
teamBin6 = teamAwinSeries.count(2)/n
teamBin7 = teamAwinSeries.count(3) / n
teamAin4 = teamBwinSeries.count(0)/n
teamAin5  = teamBwinSeries.count(1)/n
teamAin6 = teamBwinSeries.count(2)/n
teamAin7 = teamBwinSeries.count(3) / n

Answer 1

This can be done easily with numpy (Python 2.7)

import numpy as np

probs = np.array([.5 ,.56 ,.47 ,.55 ,.59 ,.59 ,.38])
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()}

Output:

>>> sum(teamAwinprops.values()) + sum(teamBwinprops.values())
1.0
>>> teamAwincount
{4: 26186, 5: 47062, 6: 59222, 7: 95381}
>>> teamBwincount
{4: 36187, 5: 79695, 6: 97802, 7: 58465}