Code to answer https://rpg.stackexchange.com/q/146809/40516

statroll.R

library(dplyr)
library(matrixStats)
library(stringr)
library(ggplot2)

## Using point buy table from https://www.d20pfsrd.com/basics-ability-scores/ability-scores/
## Assume anything less than 7 is just -4 points
point_buy_values <-
    c(-4, -4, -4, -4, -4, -4, -4, -2, -1, 0, 1, 2, 3, 5, 7, 10, 13, 17) %>%
    set_names(seq_along(.))

print(point_buy_values)

d6 <- function(n=1) sample.int(n = 6, size = n, replace = TRUE)

# 4d6, drop lowest
roll_stat <- function(n=1) {
    all_rolls <- matrix(d6(n*4), ncol=4)
    rowSums(all_rolls) - rowMins(all_rolls)
}

# Function to roll a 6x6 matrix of stats
roll_stat_matrix <- function() {
    matrix(roll_stat(n = 6*6), nrow=6) %>%
        set_rownames(str_c("R", seq_len(nrow(.)))) %>%
        set_colnames(str_c("C", seq_len(ncol(.))))
}

# Function to get the rows, columns, and diagonals
get_stat_arrays <- function(statmat) {
    rbind(statmat, t(statmat),
          LR=diag(statmat),
          RL=statmat %>% .[row(.) + col(.) - nrow(.) == 1]) %>%
        set_colnames(NULL)
}

# Function to convert stats to point buy values
get_pb <- function(x) {
    y <- x
    y[] <- point_buy_values[x]
    y
}

# Function to check the rows, columns, and diagonals for the best point buy total
get_best_pb_from_stat_matrix <- function(x) {
    x %>%
        get_stat_arrays() %>%
        get_pb() %>%
        rowSums() %>%
        max
}

# Let's generate a bunch of stat matrices and find the best point buy value for each one
best_pb <- replicate(n=500000, expr=get_best_pb_from_stat_matrix(roll_stat_matrix()))

# Now we print out the mean, median, etc., and plot the histogram
ggplot(data_frame(best_pb = best_pb)) + 
    aes(x=best_pb) + 
    geom_histogram(aes(y=..count../sum(..count..)), 
                   binwidth = 1, center = 0) + 
    xlab("Best point-buy value") + 
    ylab("Frequency")
print(summary(best_pb))