API tools faq
paste
justsing

MPQ 4 star pulling v2

justsing
Jun 6th, 2018
131
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
R 19.50 KB | None | 0 0
raw download clone embed print report
  1. ### title: MPQ 4* pulling v2
  2. ### author: justsing
  3. ### date: June 7, 2018
  4.  
  5. ## Initial Parameters
  6.  
  7. library("data.table")
  8. library("ggplot2")
  9. library("plyr")
  10.  
  11. set.seed(17)
  12. nsims = 10000
  13. options(digits = 3)
  14.  
  15. n4 = 65
  16. nlatest = 12
  17. n5 = 3
  18. p4 = 0.85
  19. p_old_latest = p4/(2*nlatest)
  20. p_old_vintage = p4/(2*(n4-nlatest))
  21. p_new = p4/n4
  22. p_BH = 0.05
  23.  
  24. plist_old = c(rep(p_old_latest, nlatest), rep(p_old_vintage, (n4-nlatest)))
  25. plist_new = rep(p_new, n4)
  26. plist_five = rep((1-p4)/n5, n5)
  27.  
  28.  
  29. ## Function for finding number of pulls needed
  30.  
  31. pull_lt <- function(plist4, plist5, specific, BH, old, chosen){
  32.   totalpulls = 0
  33.   fours = matrix(0, nrow=n4, ncol=3)
  34.   four_numcovers = rep(0, n4)
  35.   four_maxcovers = 0
  36.   while (four_maxcovers < 13){
  37.     totalpulls = totalpulls + 1
  38.     pulled_four = rbinom(1, 1, p4)
  39.     pulled_BH = rbinom(1, 1, p_BH)
  40.     if (pulled_four >= 1){
  41.       char = sample(1:n4, 1, prob=plist4)
  42.       cover = sample(1:3, 1, prob=c(1, 1, 1))
  43.       fours[char, cover] = min((fours[char, cover]+1*(four_numcovers[char] < 13)), 5)
  44.       four_numcovers = rowSums(fours)
  45.       if (BH >= 1 & pulled_BH >= 1){
  46.         cover = sample(1:3, 1, prob=c(1, 1, 1))
  47.         fours[chosen, cover] = min((fours[chosen, cover]+1*(four_numcovers[chosen] < 13)), 5)
  48.         four_numcovers = rowSums(fours)
  49.       }
  50.     }
  51.     four_maxcovers = ifelse(specific >= 1, four_numcovers[chosen], max(four_numcovers))
  52.   }
  53.   four_dist = rep(0, 14)
  54.   four_dist[2:14] = tabulate(four_numcovers, 13)
  55.   four_dist[1] = n4 - sum(four_dist[2:14])
  56.   c(old, specific, BH, chosen, totalpulls, four_dist)
  57. }
  58.  
  59.  
  60.  
  61. ## Run simulations for pulls needed
  62.  
  63. ptm <- proc.time()
  64.  
  65. old_latest_noBH = t(replicate(nsims, pull_lt(plist4=plist_old, plist5=plist_five, specific=1, BH=0, old=1, chosen=1)))
  66.  
  67. old_latest_BH = t(replicate(nsims, pull_lt(plist4=plist_old, plist5=plist_five, specific=1, BH=1, old=1, chosen=1)))
  68.  
  69. old_vintage_noBH = t(replicate(nsims, pull_lt(plist4=plist_old, plist5=plist_five, specific=1, BH=0, old=1, chosen=n4)))
  70.  
  71. old_vintage_BH = t(replicate(nsims, pull_lt(plist4=plist_old, plist5=plist_five, specific=1, BH=1, old=1, chosen=n4)))
  72.  
  73. old_any_noBH = t(replicate(nsims, pull_lt(plist4=plist_old, plist5=plist_five, specific=0, BH=0, old=1, chosen=1)))
  74.  
  75. old_any_BH = t(replicate(nsims, pull_lt(plist4=plist_old, plist5=plist_five, specific=0, BH=1, old=1, chosen=1)))
  76.  
  77. old_any_BH_vint = t(replicate(nsims, pull_lt(plist4=plist_old, plist5=plist_five, specific=0, BH=1, old=1, chosen=n4)))
  78.  
  79.  
  80. new_latest_noBH = t(replicate(nsims, pull_lt(plist4=plist_new, plist5=plist_five, specific=1, BH=0, old=0, chosen=1)))
  81.  
  82. new_latest_BH = t(replicate(nsims, pull_lt(plist4=plist_new, plist5=plist_five, specific=1, BH=1, old=0, chosen=1)))
  83.  
  84. new_vintage_noBH = t(replicate(nsims, pull_lt(plist4=plist_new, plist5=plist_five, specific=1, BH=0, old=0, chosen=n4)))
  85.  
  86. new_vintage_BH = t(replicate(nsims, pull_lt(plist4=plist_new, plist5=plist_five, specific=1, BH=1, old=0, chosen=n4)))
  87.  
  88. new_any_noBH = t(replicate(nsims, pull_lt(plist4=plist_new, plist5=plist_five, specific=0, BH=0, old=0, chosen=1)))
  89.  
  90. new_any_BH = t(replicate(nsims, pull_lt(plist4=plist_new, plist5=plist_five, specific=0, BH=1, old=0, chosen=1)))
  91.  
  92. new_any_BH_vint = t(replicate(nsims, pull_lt(plist4=plist_new, plist5=plist_five, specific=0, BH=1, old=0, chosen=n4)))
  93.  
  94.  
  95. pull_results = as.data.frame(rbind(old_latest_noBH, old_latest_BH, old_vintage_noBH, old_vintage_BH, old_any_noBH, old_any_BH, old_any_BH_vint, new_latest_noBH, new_latest_BH, new_vintage_noBH, new_vintage_BH, new_any_noBH, new_any_BH, new_any_BH_vint))
  96.  
  97. colnames(pull_results) = c("old", "specific", "BH", "chosen", "totalpulls", "covers00", "covers01", "covers02", "covers03", "covers04", "covers05", "covers06", "covers07", "covers08", "covers09", "covers10", "covers11", "covers12", "covers13")
  98.  
  99. proc.time() - ptm
  100.  
  101.  
  102. ## Results for pulls needed
  103.  
  104. # Results for Total Pulls Needed
  105. melt(setDT(pull_results[, c("old", "specific", "BH", "chosen", "totalpulls")]), id=c("specific", "chosen", "old", "BH"))[, .(mean = mean(value),
  106.                                                              sd = sd(value),
  107.                                                              median = median(value),
  108.                                                              p25 = quantile(value, probs=0.25),
  109.                                                              p75 = quantile(value, probs=0.75),
  110.                                                              p5 = quantile(value, probs=0.05),
  111.                                                              p95 = quantile(value, probs=0.95),
  112.                                                              min = min(value),
  113.                                                              max = max(value)),
  114.                                                          .(specific, chosen, old, BH, variable)]
  115.  
  116. # Boxplot of Total Pulls Needed
  117.  
  118. colorgroups = factor(pull_results$specific*pull_results$chosen,
  119. levels = c(0, 1, n4),
  120. labels = c("cover any 4*", "cover specific latest", "cover specific vintage"))
  121.  
  122. ggplot(pull_results, aes(x=interaction(specific, chosen, old, BH), y=totalpulls, color=colorgroups)) +
  123.   geom_boxplot() +
  124.   coord_trans(y = "log10") +
  125.   scale_x_discrete(limits=c("0.1.1.0", "0.1.0.0", "0.1.1.1", "0.1.0.1", "0.65.1.1", "0.65.0.1", "1.1.1.0", "1.1.0.0", "1.1.1.1", "1.1.0.1", "1.65.1.0", "1.65.0.0", "1.65.1.1", "1.65.0.1"), labels=c("0.1.1.0"="old", "0.1.0.0"="new", "0.1.1.1"="old + latest BH", "0.1.0.1"="new + latest BH", "0.65.1.1"="old + vintage BH", "0.65.0.1"="new + vintage BH", "1.1.1.0" = "old", "1.1.0.0"="new", "1.1.1.1" = "old + latest BH", "1.1.0.1" = "new + latest BH", "1.65.1.0" = "old", "1.65.0.0" = "new", "1.65.1.1" = "old + vintage BH", "1.65.0.1" = "new + vintage BH")) +
  126.   xlab("4* Odds in LTs") +
  127.   ylab("Number of Pulls Needed to Cover a 4*") +
  128.   ggtitle("Number of Pulls Needed to Cover Any 4* vs. a Specific Latest or Vintage 4*") +
  129.   theme(legend.position=c(0.15, 0.82), legend.title=element_blank(), legend.background = element_rect(fill=alpha('white', 0.6)), axis.text.x=element_text(angle=50, hjust=1, size=rel(1.1)))
  130.  
  131. # Results on 4* coverage
  132.  
  133. setDT(pull_results)[, as.list(unlist(lapply(.SD, function(x) list(mean = mean(x),
  134.                                                             sd = sd(x))))),
  135.               by = .(specific, chosen, old, BH), .SDcols=c("covers13","covers12")]
  136.  
  137.  
  138. ## Function for covered 4*s based on set number of pulls
  139.  
  140. set_pull_lt <- function(plist4, plist5, BH, old, chosen, num_pulls){
  141.   chosen0 = chosen
  142.   fours = matrix(0, nrow=n4, ncol=3)
  143.   four_numcovers = rep(0, n4)
  144.   for (totalpulls in 1:num_pulls){
  145.     pulled_four = rbinom(1, 1, p4)
  146.     pulled_BH = rbinom(1, 1, p_BH)
  147.     if (pulled_four >= 1){
  148.       char = sample(1:n4, 1, prob=plist4)
  149.       cover = sample(1:3, 1, prob=c(1, 1, 1))
  150.       fours[char, cover] = min((fours[char, cover]+1*(four_numcovers[char] < 13)), 5)
  151.       four_numcovers = rowSums(fours)
  152.       if (BH >= 1 & pulled_BH >= 1){
  153.         cover = sample(1:3, 1, prob=c(1, 1, 1))
  154.         fours[chosen, cover] = min((fours[chosen, cover]+1*(four_numcovers[chosen] < 13)), 5)
  155.         four_numcovers = rowSums(fours)
  156.       }
  157.     }
  158.    
  159.     if (four_numcovers[chosen] >= 13){chosen = chosen + 1}
  160.   }
  161.   four_dist = rep(0, 14)
  162.   four_dist[2:14] = tabulate(four_numcovers, 13)
  163.   four_dist[1] = n4 - sum(four_dist[2:14])
  164.  
  165.   num_covers_specific = c(four_numcovers[1], four_numcovers[n4])
  166.   covered_latest = sum(four_numcovers[1:nlatest] >= 13)
  167.   covered_vintage = sum(four_numcovers[(nlatest+1):n4] >= 13)
  168.   c(old, BH, chosen0, totalpulls, num_covers_specific, covered_latest, covered_vintage, four_dist)
  169. }
  170.  
  171.  
  172. ## Run simulations for set number of pulls
  173.  
  174. ptm <- proc.time()
  175.  
  176. set_old_noBH_300 = t(replicate(nsims, set_pull_lt(plist4=plist_old, plist5=plist_five, BH=0, old=1, chosen=1, num_pulls=300)))
  177.  
  178. set_old_noBH_600 = t(replicate(nsims, set_pull_lt(plist4=plist_old, plist5=plist_five, BH=0, old=1, chosen=1, num_pulls=600)))
  179.  
  180. set_old_latestBH_300 = t(replicate(nsims, set_pull_lt(plist4=plist_old, plist5=plist_five, BH=1, old=1, chosen=1, num_pulls=300)))
  181.  
  182. set_old_vintageBH_300 = t(replicate(nsims, set_pull_lt(plist4=plist_old, plist5=plist_five, BH=1, old=1, chosen=(nlatest+1), num_pulls=300)))
  183.  
  184. set_old_latestBH_600 = t(replicate(nsims, set_pull_lt(plist4=plist_old, plist5=plist_five, BH=1, old=1, chosen=1, num_pulls=600)))
  185.  
  186. set_old_vintageBH_600 = t(replicate(nsims, set_pull_lt(plist4=plist_old, plist5=plist_five, BH=1, old=1, chosen=(nlatest+1), num_pulls=600)))
  187.  
  188.  
  189. set_new_noBH_300 = t(replicate(nsims, set_pull_lt(plist4=plist_new, plist5=plist_five, BH=0, old=0, chosen=1, num_pulls=300)))
  190.  
  191. set_new_noBH_600 = t(replicate(nsims, set_pull_lt(plist4=plist_new, plist5=plist_five, BH=0, old=0, chosen=1, num_pulls=600)))
  192.  
  193. set_new_latestBH_300 = t(replicate(nsims, set_pull_lt(plist4=plist_new, plist5=plist_five, BH=1, old=0, chosen=1, num_pulls=300)))
  194.  
  195. set_new_vintageBH_300 = t(replicate(nsims, set_pull_lt(plist4=plist_new, plist5=plist_five, BH=1, old=0, chosen=(nlatest+1), num_pulls=300)))
  196.  
  197. set_new_latestBH_600 = t(replicate(nsims, set_pull_lt(plist4=plist_new, plist5=plist_five, BH=1, old=0, chosen=1, num_pulls=600)))
  198.  
  199. set_new_vintageBH_600 = t(replicate(nsims, set_pull_lt(plist4=plist_new, plist5=plist_five, BH=1, old=0, chosen=(nlatest+1), num_pulls=600)))
  200.  
  201.  
  202. set_pull_results = as.data.frame(rbind(set_old_noBH_300, set_old_latestBH_300, set_old_vintageBH_300, set_new_noBH_300, set_new_latestBH_300, set_new_vintageBH_300, set_old_noBH_600, set_old_latestBH_600, set_old_vintageBH_600, set_new_noBH_600, set_new_latestBH_600, set_new_vintageBH_600))
  203.  
  204. colnames(set_pull_results) = c("old", "BH", "chosen", "totalpulls", "num_covers_specific_latest", "num_covers_specific_vintage", "covered_latest", "covered_vintage", "covers00", "covers01", "covers02", "covers03", "covers04", "covers05", "covers06", "covers07", "covers08", "covers09", "covers10", "covers11", "covers12", "covers13")
  205.  
  206. set_pull_results$covered_specific_latest = 1*(set_pull_results$num_covers_specific_latest >= 13)
  207.  
  208. set_pull_results$covered_specific_vintage = 1*(set_pull_results$num_covers_specific_vintage >= 13)
  209.  
  210. set_pull_results$covers_ge10 = set_pull_results$covers10 + set_pull_results$covers11 + set_pull_results$covers12 + set_pull_results$covers13
  211.  
  212. proc.time() - ptm
  213.  
  214.  
  215. ## Results for set pulls
  216.  
  217. # Results on 4* coverage
  218.  
  219. # Number of Specific Latest Covers
  220. melt(setDT(set_pull_results[, c("old", "BH", "chosen", "totalpulls", "num_covers_specific_latest")]), id=c("totalpulls", "BH", "old", "chosen"))[, .(mean = mean(value),
  221.                                                              sd = sd(value),
  222.                                                              median = median(value),
  223.                                                              p25 = quantile(value, probs=0.25),
  224.                                                              p75 = quantile(value, probs=0.75),
  225.                                                              p5 = quantile(value, probs=0.05),
  226.                                                              p95 = quantile(value, probs=0.95),
  227.                                                              min = min(value),
  228.                                                              max = max(value)),
  229.                                                          .(totalpulls, old, BH, chosen, variable)]
  230.  
  231. # Number of Specific Vintage Covers
  232. melt(setDT(set_pull_results[, c("old", "BH", "chosen", "totalpulls", "num_covers_specific_vintage")]), id=c("totalpulls", "BH", "old", "chosen"))[, .(mean = mean(value),
  233.                                                              sd = sd(value),
  234.                                                              median = median(value),
  235.                                                              p25 = quantile(value, probs=0.25),
  236.                                                              p75 = quantile(value, probs=0.75),
  237.                                                              p5 = quantile(value, probs=0.05),
  238.                                                              p95 = quantile(value, probs=0.95),
  239.                                                              min = min(value),
  240.                                                              max = max(value)),
  241.                                                          .(totalpulls, old, BH, chosen, variable)]
  242.  
  243. # Probability of covering specific latest or vintage
  244. setDT(set_pull_results)[, as.list(unlist(lapply(.SD, function(x) list(mean = mean(x),
  245.                                                             sd = sd(x)
  246.                                                             )))),
  247.               by = .(totalpulls, old, BH, chosen), .SDcols=c("covered_specific_latest", "covered_specific_vintage")]
  248.  
  249. # Number of covered 4*s
  250.  
  251. melt(setDT(set_pull_results[, c("old", "BH", "chosen", "totalpulls", "covers13")]), id=c("totalpulls", "BH", "old", "chosen"))[, .(mean = mean(value),
  252.                                                              sd = sd(value),
  253.                                                              median = median(value),
  254.                                                              p25 = quantile(value, probs=0.25),
  255.                                                              p75 = quantile(value, probs=0.75),
  256.                                                              p5 = quantile(value, probs=0.05),
  257.                                                              p95 = quantile(value, probs=0.95),
  258.                                                              min = min(value),
  259.                                                              max = max(value)),
  260.                                                          .(totalpulls, old, BH, chosen, variable)]
  261.  
  262. melt(setDT(set_pull_results[, c("old", "BH", "chosen", "totalpulls", "covered_latest")]), id=c("totalpulls", "BH", "old", "chosen"))[, .(mean = mean(value),
  263.                                                              sd = sd(value),
  264.                                                              median = median(value),
  265.                                                              p25 = quantile(value, probs=0.25),
  266.                                                              p75 = quantile(value, probs=0.75),
  267.                                                              p5 = quantile(value, probs=0.05),
  268.                                                              p95 = quantile(value, probs=0.95),
  269.                                                              min = min(value),
  270.                                                              max = max(value)),
  271.                                                          .(totalpulls, old, BH, chosen, variable)]
  272.  
  273. melt(setDT(set_pull_results[, c("old", "BH", "chosen", "totalpulls", "covered_vintage")]), id=c("totalpulls", "BH", "old", "chosen"))[, .(mean = mean(value),
  274.                                                              sd = sd(value),
  275.                                                              median = median(value),
  276.                                                              p25 = quantile(value, probs=0.25),
  277.                                                              p75 = quantile(value, probs=0.75),
  278.                                                              p5 = quantile(value, probs=0.05),
  279.                                                              p95 = quantile(value, probs=0.95),
  280.                                                              min = min(value),
  281.                                                              max = max(value)),
  282.                                                          .(totalpulls, old, BH, chosen, variable)]
  283.  
  284. # Number of 4*s with 10 or more covers
  285.  
  286. melt(setDT(set_pull_results[, c("old", "BH", "chosen", "totalpulls", "covers_ge10")]), id=c("totalpulls", "BH", "old", "chosen"))[, .(mean = mean(value),
  287.                                                              sd = sd(value),
  288.                                                              median = median(value),
  289.                                                              p25 = quantile(value, probs=0.25),
  290.                                                              p75 = quantile(value, probs=0.75),
  291.                                                              p5 = quantile(value, probs=0.05),
  292.                                                              p95 = quantile(value, probs=0.95),
  293.                                                              min = min(value),
  294.                                                              max = max(value)),
  295.                                                          .(totalpulls, old, BH, chosen, variable)]
  296.  
  297. # Box Plot of covered 4*s
  298.  
  299. ggplot(set_pull_results, aes(x=interaction(old, BH, chosen), y=covers13, color=as.factor(interaction(old, BH, chosen)))) +
  300.   geom_boxplot() +
  301.   facet_wrap(~paste(totalpulls, "pulls")) +
  302.   scale_x_discrete(limits=c("1.0.1", "0.0.1", "1.1.1", "0.1.1", "1.1.13", "0.1.13"),labels=c("1.0.1" = "old", "0.0.1" = "new", "1.1.1" = "old + latest BH", "0.1.1" = "new + latest BH", "1.1.13" = "old + vintage BH", "0.1.13" = "new + vintage BH")) +
  303.   xlab("4* Odds in LTs") +
  304.   ylab("Number of Covered 4*s") +
  305.   ggtitle("Number of Covered 4*s in a Set Number of Pulls") +
  306.   theme(legend.position="none", legend.title=element_blank(), axis.text.x=element_text(angle=50, hjust=1, size=rel(1.1)))
  307.  
  308. # Box Plot of 4*s with at least 10 covers
  309.  
  310. ggplot(set_pull_results, aes(x=interaction(old, BH, chosen), y=covers_ge10, color=as.factor(interaction(old, BH, chosen)))) +
  311.   geom_boxplot() +
  312.   facet_wrap(~paste(totalpulls, "pulls")) +
  313.   scale_x_discrete(limits=c("1.0.1", "0.0.1", "1.1.1", "0.1.1", "1.1.13", "0.1.13"),labels=c("1.0.1" = "old", "0.0.1" = "new", "1.1.1" = "old + latest BH", "0.1.1" = "new + latest BH", "1.1.13" = "old + vintage BH", "0.1.13" = "new + vintage BH")) +
  314.   xlab("4* Odds in LTs") +
  315.   ylab("Number of 4*s with 10 or More Covers") +
  316.   ggtitle("Number of 4*s with 10 or More Covers in a Set Number of Pulls") +
  317.   theme(legend.position="none", legend.title=element_blank(), axis.text.x=element_text(angle=50, hjust=1, size=rel(1.1)))
  318.  
  319. # Bar Chart of Mean Covered 4*s by Type
  320.  
  321. data_summary <- function(data, varname, groupnames){
  322.   require(plyr)
  323.   summary_func <- function(x, col){
  324.     c(varmean = mean(x[[col]], na.rm=TRUE),
  325.       sd = sd(x[[col]], na.rm=TRUE))
  326.   }
  327.   data_sum<-ddply(data, groupnames, .fun=summary_func,
  328.                   varname)
  329.   #data_sum <- rename(data_sum, c("mean" = varname))
  330.  return(data_sum)
  331. }
  332.  
  333. set_pull_summary_latest <- data_summary(set_pull_results, varname = "covered_latest", groupnames=c("totalpulls", "old", "BH", "chosen"))
  334.  
  335. set_pull_summary_latest$fourtype = "latest"
  336.  
  337. set_pull_summary_vintage <- data_summary(set_pull_results, varname = "covered_vintage", groupnames=c("totalpulls", "old", "BH", "chosen"))
  338.  
  339. set_pull_summary_vintage$fourtype = "vintage"
  340.  
  341. set_pull_summary = rbind(set_pull_summary_latest, set_pull_summary_vintage)
  342.  
  343. ggplot(set_pull_summary, aes(x=interaction(old, BH, chosen), y=varmean, fill=fourtype)) +
  344.    geom_bar(stat="identity") +
  345.   facet_wrap(~paste(totalpulls, "pulls")) +
  346.   scale_x_discrete(limits=c("1.0.1", "0.0.1", "1.1.1", "0.1.1", "1.1.13", "0.1.13"),labels=c("1.0.1" = "old", "0.0.1" = "new", "1.1.1" = "old + latest BH", "0.1.1" = "new + latest BH", "1.1.13" = "old + vintage BH", "0.1.13" = "new + vintage BH")) +
  347.   xlab("Odds for 4* in LT") +
  348.   ylab("Mean Number of Covered 4*s") +
  349.   ggtitle("Mean Number of Covered 4*s by Type in a Set Number of Pulls") +
  350.   theme(legend.title=element_blank(), axis.text.x=element_text(angle=50, hjust=1, size=rel(1.1)))
Advertisement
Add Comment
Please, Sign In to add comment
Public Pastes
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!