heroprob <- 0.14

prob <- c(0.84, 0.15, 0.01)
values <- c(64, 256, 512)

ncalls <- 10^7

card <- function(n, heroprob, prob, values) {
  ifelse(runif(n) <= heroprob,
         sample(values, n, replace = TRUE, prob=prob),
         0)
}

call <- function(n, heroprob, prob, values) {
  cbind(card(n, heroprob, prob, values),
        card(n, heroprob, prob, values),
        card(n, heroprob, prob, values))
}


calls <- call(ncalls, heroprob, prob, values)
callvalues <- pmax(calls[,1], calls[,2], calls[,3])
sum(callvalues) / ncalls