Script to answer https://rpg.stackexchange.com/questions/156379

colville-rollstat.R

## License: WTFPL
library(reshape2)
library(matrixStats)
library(ggplot2)
library(magrittr)
library(dplyr)
## You also need the Hmisc package installed, but I don't load it
## because of namespace clashes with dplyr

## Enumerate all possible 4d6 rolls, then order each row from high to
## low. Each row has an equal probability of being rolled
all_stat_rolls <- expand.grid(1:6, 1:6, 1:6, 1:6) %>%
    as.matrix() %>%
    rowQuantiles(probs = seq(1, 0, length.out=4)) %>%
    set_colnames(paste0("d6_", LETTERS[1:4])) %>%
    as_tibble %>%
    mutate(Stat = d6_A + d6_B + d6_C)

all_stat_rolls_min8 <- all_stat_rolls %>% filter(Stat >= 8)

## Probability of rolling each stat from 8 to 18.
min8_stat_dist <- all_stat_rolls_min8 %>%
    group_by(Stat) %>% dplyr::summarize(Prob = length(Stat) / nrow(all_stat_rolls_min8))

## Return probability of rolling x via 4d6 keep 3 highest (it just
## looks up the values in the above table)
prob_4d6kh3 <- function(x) {
    min8_stat_dist %$% {
        Prob[match(x, Stat)]
    }
}

## Every possible stat roll, with the probability of rolling it
all_stat_sets <- do.call(expand.grid, rep(list(8:18), 6)) %>%
    as.matrix %>%
    rowQuantiles(probs=seq(1,0,length.out=6)) %>%
    round %>% as_tibble %>%
    set_colnames(LETTERS[1:6]) %>%
    mutate(Prob = Reduce(`*`, Map(prob_4d6kh3, .[LETTERS[1:6]]))) %>%
    group_by(A, B, C, D, E, F) %>%
    dplyr::summarize(Prob = sum(Prob)) %>%
    ungroup
## Total probability should be 1 (plus or minus floating point rounding error)
stopifnot(abs(sum(all_stat_sets$Prob) - 1) <= .Machine$double.eps * 100)

## Select only Colville-approved stat sets (at least 2 stats 15 or
## higher) and renormalize the probabilities to sum to 1
all_colville_stat_sets <-
    all_stat_sets %>%
    filter(B >= 15) %>%
    mutate(Prob = Prob / sum(Prob))
stopifnot(abs(sum(all_colville_stat_sets$Prob) - 1) <= .Machine$double.eps * 100)

## Extract the probability distribution for each stat
all_stat_probs <- do.call(rbind, lapply(LETTERS[1:6], function(i) {
    all_colville_stat_sets %>%
        rename(Value = !!i) %>%
        group_by(Value) %>%
        dplyr::summarize(Prob=sum(Prob)) %>%
        cbind(Stat = i, .)
})) %>%
    ## Fill in 0 probabilities for stat rolls that never occur
    ## (because they're filtered out), for plotting purposes
    right_join(expand.grid(Stat=LETTERS[1:6], Value=8:18)) %>%
    mutate(Prob = ifelse(is.na(Prob), 0, Prob))

## Make a nice-looking plot of the distribution of each stat
p <- ggplot(all_stat_probs) +
    geom_ribbon(aes(x = Value, ymax = Prob*100, fill = Stat),
                color = NA, ymin=0, alpha=0.1) +
    geom_line(aes(x = Value, y = Prob*100, color = Stat)) +
    scale_x_continuous(limits = range(min8_stat_dist$Stat),
                       breaks = sort(min8_stat_dist$Stat),
                       minor_breaks = NULL) +
    ggtitle("Colville stat distributions",
            "A is best stat, B is 2nd best, etc.") +
    xlab("Stat Value") +
    ylab("Percent chance")
print(p)

## Print out mean, median, etc. for each of the 6 stats
all_stat_probs %>%
    filter(Prob > 0) %>%
    group_by(Stat) %>%
    dplyr::summarize(
        Min = min(Value),
        Q25 = Hmisc::wtd.quantile(Value, Prob, .25, normwt = TRUE),
        Median = Hmisc::wtd.quantile(Value, Prob, .5, normwt = TRUE),
        Mean = Hmisc::wtd.mean(Value, Prob, normwt = TRUE),
        Q75 = Hmisc::wtd.quantile(Value, Prob, .75, normwt = TRUE),
        Max = max(Value))