Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
## 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))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.