MPQ 4 star pulling

### MPQ 4* pulling
### justsing
### June 6, 2018

## Setup

n4 = 65
nlatest = 12
p4_lt = 0.85
p_lt_current = p4_lt/(2*nlatest*3)
p_lt_new = p4_lt/(n4*3)

coptions = c(1, 2, 3, 4, 5)

## Simulation Functions

pull_noBH <- function(coverlist, pcover, p4){
  totalpulls = 0
  covers = rep(0, 5)
  plist = c(rep(pcover, 3), (p4-3*pcover), (1-p4))
  while (sum(covers[1:3]) < 13){
    totalpulls = totalpulls + 1
    pull = sample(coverlist, 1, prob=plist)
    covers[pull] = min((covers[pull]+1), 5)
  }
  totalpulls
}

pull_BH <- function(coverlist, pcover, p4){
  totalpulls = 0
  covers = rep(0, 5)
  plist = c(rep(pcover, 3), (p4-3*pcover), (1-p4))
  while (sum(covers[1:3]) < 13){
    totalpulls = totalpulls + 1
    pull = sample(coverlist, 1, prob=plist)
    covers[pull] = min((covers[pull]+1), 5)
    if (pull <= 4){
      bonus = sample(coverlist[1:4], 1, prob=c(0.05/3, 0.05/3, 0.05/3, 0.95))
      covers[bonus] = min((covers[bonus]+1), 5)
    }
  }
  totalpulls
}

### Run simulations

set.seed(17)
nsims = 100000

lt_current_noBH = replicate(nsims, pull_noBH(coverlist=coptions, pcover=p_lt_current, p4=p4_lt))
lt_current_BH = replicate(nsims, pull_BH(coverlist=coptions, pcover=p_lt_current, p4=p4_lt))
lt_new_noBH = replicate(nsims, pull_noBH(coverlist=coptions, pcover=p_lt_new, p4=p4_lt))
lt_new_BH = replicate(nsims, pull_BH(coverlist=coptions, pcover=p_lt_new, p4=p4_lt))

lt_results = as.data.frame(cbind(lt_current_noBH, lt_current_BH, lt_new_noBH, lt_new_BH))
summary(lt_results)

sapply(lt_results, mean)
sapply(lt_results, sd)
sapply(lt_results, quantile, probs = c(0.5, 0.25, 0.75, 0.05, 0.95))