# assuming known team distributions, this code # simulates the performance of two competing units # number of sims to run nsims <- 10000 # possessions in each game nposs <- 100 # beta distribution class setClass("beta_params", representation(alpha="numeric",beta="numeric")); # unit class setClass("unit", representation( # offensive rates # distribution of 2pt shots, 3pt shots, and turnovers on each play o_2pt = "numeric", o_3pt = "numeric", o_tov = "numeric", # FG% for 2pt and 3pt shots o_fgp_2pt = "beta_params", o_fgp_3pt = "beta_params", # distribution of 2pt and 3pt shots that have a shooting foul o_sf_2pt = "beta_params", o_sf_3pt = "beta_params", # free throw % o_ftp = "beta_params", # offensive rebound % o_rebp = "beta_params", # defensive rates # distribution of 2pt shots, 3pt shots, and turnovers on each play d_2pt = "numeric", d_3pt = "numeric", d_tov = "numeric", # FG% for 2pt and 3pt shots d_fgp_2pt = "beta_params", d_fgp_3pt = "beta_params", # distribution of 2pt and 3pt shots that have a shooting foul d_sf_2pt = "beta_params", d_sf_3pt = "beta_params", # defensive rebound % d_rebp = "beta_params" )) # add in prior information # for this case, we assume a uniform prior, which simply means we add 1 to all params "add_prior" <- function(unit) { unit@o_2pt <- 1+unit@o_2pt unit@o_3pt <- 1+unit@o_3pt unit@o_tov <- 1+unit@o_tov unit@o_fgp_2pt@alpha <- 1+unit@o_fgp_2pt@alpha unit@o_fgp_2pt@beta <- 1+unit@o_fgp_2pt@beta unit@o_fgp_3pt@alpha <- 1+unit@o_fgp_3pt@alpha unit@o_fgp_3pt@beta <- 1+unit@o_fgp_3pt@beta unit@o_sf_2pt@alpha <- 1+unit@o_sf_2pt@alpha unit@o_sf_2pt@beta <- 1+unit@o_sf_2pt@beta unit@o_sf_3pt@alpha <- 1+unit@o_sf_3pt@alpha unit@o_sf_3pt@beta <- 1+unit@o_sf_3pt@beta unit@o_ftp@alpha <- 1+unit@o_ftp@alpha unit@o_ftp@beta <- 1+unit@o_ftp@beta unit@o_rebp@alpha <- 1+unit@o_rebp@alpha unit@o_rebp@beta <- 1+unit@o_rebp@beta unit@d_2pt <- 1+unit@d_2pt unit@d_3pt <- 1+unit@d_3pt unit@d_tov <- 1+unit@d_tov unit@d_fgp_2pt@alpha <- 1+unit@d_fgp_2pt@alpha unit@d_fgp_2pt@beta <- 1+unit@d_fgp_2pt@beta unit@d_fgp_3pt@alpha <- 1+unit@d_fgp_3pt@alpha unit@d_fgp_3pt@beta <- 1+unit@d_fgp_3pt@beta unit@d_sf_2pt@alpha <- 1+unit@d_sf_2pt@alpha unit@d_sf_2pt@beta <- 1+unit@d_sf_2pt@beta unit@d_sf_3pt@alpha <- 1+unit@d_sf_3pt@alpha unit@d_sf_3pt@beta <- 1+unit@d_sf_3pt@beta unit@d_rebp@alpha <- 1+unit@d_rebp@alpha unit@d_rebp@beta <- 1+unit@d_rebp@beta return(unit) } # draws from the dirichlet distribution # code taken from MCMCpack by Kevin Rompala, # which was taken from "Greg's Miscellaneous Functions" "rdirichlet" <- function(n, alpha) { l <- length(alpha) x <- matrix(rgamma(l*n,alpha),ncol=l,byrow=TRUE) sm <- x%*%rep(1,l) return(x/as.vector(sm)) } # handles free throws used in simulation "handle_fts" <- function(num,pr) { points <- 0 for (i in 1:num) { result <- rbinom(1,size=1,prob=pr) if (result) { # FTM points <- points+1 } last <- result } return(c(points,last)) } "run_sims" <- function(nsims,nposs,unit1,unit2) { # vectors to keep track of the efficiency for each simulation u1.eff <- c() u2.eff <- c() # keeps track of the number of "games" each unit wins u1.wins <- 0 u2.wins <- 0 # vectors to keep track of the margins of victory for each team u1.margins <- c() u2.margins <- c() for (i in 1:nsims) { # determine who has the initial possession has_pos <- rbinom(1,size=1,prob=0.5) # vectors to keep track of the points on each possession u1.points <- c() u2.points <- c() for (j in 1:(nposs*2)) { if (!has_pos) { # unit1 has possession has <- unit1 } else { # unit2 has possession has <- unit2 } # run through game logic go <- 1 points <- 0 while (go) { # this loop continues until the opposing unit obtains possession # determine what is going to happen: 2pt shot, 3pt shot, or turnover probs <- rdirichlet(1,c(has@o_2pt,has@o_3pt,has@o_tov)) result <- rmultinom(1, size=1, prob=c(probs[1,1],probs[1,2],probs[1,3])) # some useful variables reb <- 0 # handle result if (result[1,1]) { # 2pt shot # determine if the shot was made probs <- rbeta(1,shape1=has@o_fgp_2pt@alpha,shape2=has@o_fgp_2pt@beta) result <- rbinom(1,size=1,prob=probs) if (result) { # make; add 2 points points <- points+2 reb <- 0 go <- 0 } else { # miss reb <- 1 go <- 1 } # determine if there was a foul on the shot probs <- rbeta(1,shape1=has@o_sf_2pt@alpha,shape2=has@o_sf_2pt@beta) result <- rbinom(1,size=1,prob=probs) if (result) { # handle the foul shots probs <- rbeta(1,shape1=has@o_ftp@alpha,shape2=has@o_ftp@beta) result <- handle_fts(num=2,pr=probs) points <- points+result[1] if (!result[2]) { # the last was missed; handle rebound to determine possession reb <- 1 } else { # the last was made; change of possession reb <- 0 go <- 0 } } if (reb) { # handle rebound probs <- rbeta(1,shape1=has@o_rebp@alpha,shape2=has@o_rebp@beta) result <- rbinom(1,size=1,prob=probs) if (result) { # offensive rebound go <- 1 } else { # defensive rebound; change of possession go <- 0 } } } else if (result[2,1]) { # 3pt shot # determine if the shot was made probs <- rbeta(1,shape1=has@o_fgp_3pt@alpha,shape2=has@o_fgp_3pt@beta) result <- rbinom(1,size=1,prob=probs) if (result) { # make; add 3 points points <- points+3 reb <- 0 go <- 0 } else { # miss reb <- 1 go <- 1 } # determine if there was a foul on the shot probs <- rbeta(1,shape1=has@o_sf_3pt@alpha,shape2=has@o_sf_3pt@beta) result <- rbinom(1,size=1,prob=probs) if (result) { # handle the foul shots probs <- rbeta(1,shape1=has@o_ftp@alpha,shape2=has@o_ftp@beta) result <- handle_fts(num=3,pr=probs) points <- points+result[1] if (!result[2]) { # the last was missed; handle rebound to determine possession reb <- 1 } else { # the last was made; change of possession reb <- 0 go <- 0 } } if (reb) { # handle rebound probs <- rbeta(1,shape1=has@o_rebp@alpha,shape2=has@o_rebp@beta) result <- rbinom(1,size=1,prob=probs) if (result) { # offensive rebound go <- 1 } else { # defensive rebound; change of possession go <- 0 } } } else if (result[3,1]) { # turnover; change of possession go <- 0 } } # change possession and store data if (!has_pos) { has_pos <- 1 u1.points <- c(u1.points,points) } else { has_pos <- 0 u2.points <- c(u2.points,points) } } # end nposs sim; store efficiency u1.eff <- c(u1.eff,sum(u1.points)/nposs*100) u2.eff <- c(u2.eff,sum(u2.points)/nposs*100) if (sum(u1.points) > sum(u2.points)) { u1.wins <- u1.wins+1 } else if (sum(u2.points) > sum(u1.points)) { u2.wins <- u2.wins+1 } u1.margins <- c(u1.margins,sum(u1.points)-sum(u2.points)) u2.margins <- c(u2.margins,sum(u2.points)-sum(u1.points)) } cat("Unit1 Eff: mean=",mean(u1.eff),"; sd=",sd(u1.eff),"\n",sep="") cat("Unit2 Eff: mean=",mean(u2.eff),"; sd=",sd(u2.eff),"\n",sep="") cat("Unit1 Margins: mean=",mean(u1.margins),"; sd=",sd(u1.margins),"\n",sep="") cat("Unit2 Margins: mean=",mean(u2.margins),"; sd=",sd(u2.margins),"\n",sep="") cat("Unit1 Win%: ",u1.wins/(u1.wins+u2.wins),"\n",sep="") cat("Unit2 Win%: ",u2.wins/(u1.wins+u2.wins),"\n",sep="") } # returns beta params from known mean and standard deviation "get_beta_params" <- function(mean,sd) { k <- (1-mean)/mean alpha <- ( k/sd^2 - (k+1)^2 ) / (k+1)^3 beta <- k*alpha return(c(alpha,beta)) } "combine_units" <- function(unit1,unit2) { new_unit <- new("unit") nsims <- 50000 # play distribution is modeled using the dirichlet distribution # for this case, we can draw from each case and then use those as the # values of the new dirichlet distribution s1 <- rbeta(nsims,unit1@o_2pt,unit1@o_3pt+unit1@o_tov) s2 <- rbeta(nsims,unit2@d_2pt,unit2@d_3pt+unit2@d_tov) dist <- (s1+s2)/2 params.two <- get_beta_params(mean(dist),sd(dist)) s1 <- rbeta(nsims,unit1@o_3pt,unit1@o_2pt+unit1@o_tov) s2 <- rbeta(nsims,unit2@d_3pt,unit2@d_2pt+unit2@d_tov) dist <- (s1+s2)/2 params.three <- get_beta_params(mean(dist),sd(dist)) s1 <- rbeta(nsims,unit1@o_tov,unit1@o_2pt+unit1@o_3pt) s2 <- rbeta(nsims,unit2@d_tov,unit2@d_2pt+unit2@d_3pt) dist <- (s1+s2)/2 params.tov <- get_beta_params(mean(dist),sd(dist)) # we simply take the new alpha params as the parameters to the dirichlet distribution new_unit@o_2pt = params.two[1] new_unit@o_3pt = params.three[1] new_unit@o_tov = params.tov[1] # FG% s1 <- rbeta(nsims,unit1@o_fgp_2pt@alpha,unit1@o_fgp_2pt@beta) s2 <- rbeta(nsims,unit2@d_fgp_2pt@alpha,unit2@d_fgp_2pt@beta) dist <- (s1+s2)/2 params <- get_beta_params(mean(dist),sd(dist)) new_unit@o_fgp_2pt@alpha = params[1] new_unit@o_fgp_2pt@beta = params[2] s1 <- rbeta(nsims,unit1@o_fgp_3pt@alpha,unit1@o_fgp_3pt@beta) s2 <- rbeta(nsims,unit2@d_fgp_3pt@alpha,unit2@d_fgp_3pt@beta) dist <- (s1+s2)/2 params <- get_beta_params(mean(dist),sd(dist)) new_unit@o_fgp_3pt@alpha = params[1] new_unit@o_fgp_3pt@beta = params[2] # Pr(shooting foul) s1 <- rbeta(nsims,unit1@o_sf_2pt@alpha,unit1@o_sf_2pt@beta) s2 <- rbeta(nsims,unit2@d_sf_2pt@alpha,unit2@d_sf_2pt@beta) dist <- (s1+s2)/2 params <- get_beta_params(mean(dist),sd(dist)) new_unit@o_sf_2pt@alpha = params[1] new_unit@o_sf_2pt@beta = params[2] s1 <- rbeta(nsims,unit1@o_sf_3pt@alpha,unit1@o_sf_3pt@beta) s2 <- rbeta(nsims,unit2@d_sf_3pt@alpha,unit2@d_sf_3pt@beta) dist <- (s1+s2)/2 params <- get_beta_params(mean(dist),sd(dist)) new_unit@o_sf_3pt@alpha = params[1] new_unit@o_sf_3pt@beta = params[2] # rebounding % s1 <- rbeta(nsims,unit1@o_rebp@alpha,unit1@o_rebp@beta) s2 <- 1-rbeta(nsims,unit2@d_rebp@alpha,unit2@d_rebp@beta) dist <- (s1+s2)/2 params <- get_beta_params(mean(dist),sd(dist)) new_unit@o_rebp@alpha = params[1] new_unit@o_rebp@beta = params[2] # keep others new_unit@o_ftp@alpha = unit1@o_ftp@alpha new_unit@o_ftp@beta = unit1@o_ftp@beta # we don't care about the defensive items # d_2pt, d_3pt, d_tov, d_fgp_2pt, d_fgp_3pt, d_sf_2pt, d_sf_3pt, d_rebp cat("Play dist:",new_unit@o_2pt,new_unit@o_3pt,new_unit@o_tov,"-",new_unit@o_2pt/(new_unit@o_2pt+new_unit@o_3pt+new_unit@o_tov),new_unit@o_3pt/(new_unit@o_2pt+new_unit@o_3pt+new_unit@o_tov),new_unit@o_tov/(new_unit@o_2pt+new_unit@o_3pt+new_unit@o_tov),"\n") cat("2pt FG%:",new_unit@o_fgp_2pt@alpha,new_unit@o_fgp_2pt@beta,"-",new_unit@o_fgp_2pt@alpha/(new_unit@o_fgp_2pt@alpha+new_unit@o_fgp_2pt@beta),"\n") cat("3pt FG%:",new_unit@o_fgp_3pt@alpha,new_unit@o_fgp_3pt@beta,"-",new_unit@o_fgp_3pt@alpha/(new_unit@o_fgp_3pt@alpha+new_unit@o_fgp_3pt@beta),"\n") cat("SF, 2pt:",new_unit@o_sf_2pt@alpha,new_unit@o_sf_2pt@beta,"-",new_unit@o_sf_2pt@alpha/(new_unit@o_sf_2pt@alpha+new_unit@o_sf_2pt@beta),"\n") cat("SF, 3pt:",new_unit@o_sf_3pt@alpha,new_unit@o_sf_3pt@beta,"-",new_unit@o_sf_3pt@alpha/(new_unit@o_sf_3pt@alpha+new_unit@o_sf_3pt@beta),"\n") cat("OReb%:",new_unit@o_rebp@alpha,new_unit@o_rebp@beta,"-",new_unit@o_rebp@alpha/(new_unit@o_rebp@alpha+new_unit@o_rebp@beta),"\n") return(new_unit) } # unit definitions home <- new("unit", o_2pt = 845, o_3pt = 177, o_tov = 170, o_fgp_2pt = new("beta_params", alpha=427,beta=418), o_fgp_3pt = new("beta_params", alpha=65,beta=112), o_sf_2pt = new("beta_params", alpha=116,beta=729), o_sf_3pt = new("beta_params", alpha=0,beta=177), o_ftp = new("beta_params", alpha=194,beta=58), o_rebp = new("beta_params", alpha=119,beta=450), d_2pt = 813, d_3pt = 222, d_tov = 194, d_fgp_2pt = new("beta_params", alpha=308,beta=505), d_fgp_3pt = new("beta_params", alpha=65,beta=157), d_sf_2pt = new("beta_params", alpha=109,beta=704), d_sf_3pt = new("beta_params", alpha=5,beta=217), d_rebp = new("beta_params", alpha=352,beta=149)) away <- new("unit", o_2pt = 302, o_3pt = 95, o_tov = 59, o_fgp_2pt = new("beta_params", alpha=146,beta=156), o_fgp_3pt = new("beta_params", alpha=41,beta=54), o_sf_2pt = new("beta_params", alpha=47,beta=255), o_sf_3pt = new("beta_params", alpha=1,beta=94), o_ftp = new("beta_params", alpha=89,beta=17), o_rebp = new("beta_params", alpha=51,beta=146), d_2pt = 324, d_3pt = 83, d_tov = 65, d_fgp_2pt = new("beta_params", alpha=154,beta=170), d_fgp_3pt = new("beta_params", alpha=31,beta=52), d_sf_2pt = new("beta_params", alpha=29,beta=295), d_sf_3pt = new("beta_params", alpha=0,beta=83), d_rebp = new("beta_params", alpha=132,beta=56)) # add prior to the units home <- add_prior(home) away <- add_prior(away) # use offensive/defensive data to create new distributions home.new <- combine_units(home,away) away.new <- combine_units(away,home) # run the sims system.time(run_sims(nsims,nposs,home.new,away.new))