Getting series outcome probabilities based on individual game probabilities

Question

I have the probabilities for Team X beating Team Y in a game. Each team has two separate probabilities - the probability of winning at home and the probability of winning away.

These probabilities can then be used to calculate the expected winner of a best-of-7 series between the two team. Based on whether or not a team finished higher in the standings, they'll either play a maximum of 4 games at home, or a maximum 3. The format for these games is 2-2-1-1-1, with the higher seed playing twice at home, twice on the road, once at home (if necessary), once on the road (if necessary), and once at home (if necessary).

How could I find the probability that team one wins the series (wins 4 games) in X number (ranging from 4 to 7) of games? My initial thought was to multiply binomial probabilites, like so in R:

 games.out <- function(ph, pa) {
  df <- data.frame(w = NA, x = NA, y = NA, z = NA)                                                                                      
  for(games.played in 4:7){
    home.games <- ceiling(games.played/2)
    away.games <- floor(games.played/2)
    home.wins <- 4 - away.games
    away.wins <- away.games
    loops <- games.played - 4
    win.prob <- (dbinom(home.wins, home.games, ph)*dbinom(away.wins, away.games, pa))
    home.wins <- home.wins + 1
    away.wins <- away.wins - 1
    for(i in 1:loops) {
      temp.prob <- (dbinom(home.wins, home.games, ph)*dbinom(away.wins, away.games, pa))
      home.wins <- home.wins + 1
      away.wins <- away.wins -1
      win.prob <- temp.prob + win.prob
    }
    df[1,games.played-3] <- win.prob
    colnames(df)[games.played-3] <- paste0("Win", games.played)
  }
  return(df)
}

But, when summed, this code gave values that were much too large to make sense.

> games.out(.5, .5)
     HW4     HW5      HW6       HW7
1 0.0625 0.15625 0.234375 0.2734375
> sum(games.out(.5, .5))
[1] 0.7265625

A similar question was asked here (How to calculate probability of winning best of 7 series), but that problem simply deals with one probability associated with all seven games - here, I'm looking at two different probabilities associated with potentially different numbers of games.

Any help is definitely appreciated - I've been banging my head against the wall for a couple of hours on this.

Answer 1

Essentially you want to have a formula for the probability mass function of a negative binomial, except instead of observing the outcomes of identically distributed Bernoulli trials, we now have sequences of Bernoulli trials that are not identically distributed.

Unfortunately I do not know of a nice formula. However, we can derive a solution combinatorially by enumerating all possible paths that end when four team 1 wins or team 2 wins have been observed (these are respectively coded with 1's and 0's in the following functions). I will go through the steps by commenting in-line within these two functions in R.

get_succ_prob <- function(p, q) { # modify as necessary depending on pattern of home/away games
    # probability team 1 wins at the i-th game    
    # can extend to other game win probs
    c(p, p, q, q, p, q, p) 
}

The main work is done by this function. n_trials is the number of games in the series, n_success is the number of games required to win the series, and succ_prob is the vector (generated by the previous function), such that the i-th element is the probability that team 1 wins the i-th game.

get_win_prob <- function(n_trials, n_success, succ_prob) {

        stopifnot(length(succ_prob) == n_trials)

        # enumerate 2^n_trials sequences of 0, 1.
        tmp_mat <- as.matrix(do.call(expand.grid, rep(list(0:1), n_trials)))

        # now terminate series of games at the first instance 
        # four zeroes or ones are observed
        tmp_mat <- t(apply(tmp_mat, 1, function(x) {

            threshold_to_NA <- suppressWarnings(min(min(which(cumsum(x == 0) == n_success)), min(which(cumsum(x == 1) == n_success))))
            x[seq_len(n_trials) >= (threshold_to_NA + 1)] <- NA

            return(x)
        }
        ))

        # Find unique paths to a series ending
        tmp_mat <- unique(tmp_mat)

        # calculate probability of observing series of game outcomes
        prob_row <- 
            apply(tmp_mat, 1, function(x) prod(succ_prob * x + (1 - succ_prob) * (1 - x), 
                                               na.rm = T))

        # make data frame and aggregate by team/series length
        tmp_df <- data.frame(tmp_mat)

        tmp_df$prob <- prob_row



        tmp_df$winning_team <-
            2 - (rowSums(tmp_df[, seq_len(n_trials)], na.rm = T) == n_success)

        tmp_df$trials_required <- n_trials - rowSums(is.na(tmp_df[, seq_len(n_trials)]))


        out_df <- aggregate(tmp_df$prob, by = list(tmp_df$trials_required, tmp_df$winning_team), FUN = sum)
        names(out_df) <- c('games_required_to_win', 'winning_team', 'win_prob')

        return(out_df)
    }

Test on your example:

> get_win_prob(7, 4, get_succ_prob(0.5, 0.5))
  games_required_to_win winning_team win_prob
1                     4            1  0.06250
2                     5            1  0.12500
3                     6            1  0.15625
4                     7            1  0.15625
5                     4            2  0.06250
6                     5            2  0.12500
7                     6            2  0.15625
8                     7            2  0.15625

On more general home/away win probs:

> get_win_prob(7, 4, get_succ_prob(0.6, 0.3))
  games_required_to_win winning_team win_prob
1                     4            1 0.032400
2                     5            1 0.116640
3                     6            1 0.093960
4                     7            1 0.192456
5                     4            2 0.078400
6                     5            2 0.120960
7                     6            2 0.236880
8                     7            2 0.128304