I need a 'standard array' for a D&D-like homebrew game, but anydice chokes - how to proceed?

Question

So, I'm working on a homebrew RPG that is similar to several FRPs that use 3d6 sorts of starting values for attributes/characteristics.

I know how to calculate the standard D&D (5E or other) 'standard array' (the one it ought to have been, not what they gave us).

My problem is this: Anydice chokes (time limit of 5 seconds) on anything more than 9 characteristics and I have 12.

I have noted as you move to having more attributes, the array is scattered over more values and you end up with higher top score and lower low score than you would with a smaller pool.

Here's my anydice code:

T: {2,3,3,4,5,6}

ABILITIES: 9d (3dT)  
loop P over {1..9} {
 output P @ ABILITIES named "Ability [P]"
}

The problem is I can't expand it to 12d over 1..12 and that means I can't get (as far as I know) the array I want.

Can someone either tell me how to get this out of anydice (if I can limit accuracy a bit to speed computation and fit 12 in for instance) or explain the operations one might be able to execute in Excel or Google Sheets to compute the array.

I have 32Gb of RAM in my machine and a 12 core i5. I have no problem letting the computer chew on the problem for a time. I just don't know the steps and I can't use my usual go to (anydice) because of the 5 second limit.

Thoughts?

Answer 1

A practical solution

My Icepool Python library has a more efficient pool algorithm:

from icepool import Die

num_scores = 12
t = Die([2,3,3,4,5,6])
ability = 3 @ t
pool = ability.pool(num_scores)

for i in range(num_scores):
    output(pool[i])

You can run this script online here.

Mathematical notes

I have noted as you move to having more attributes, the array is scattered over more values and you end up with higher top score and lower low score than you would with a smaller pool.

The technical term for this is order statistics. For a continuous uniform distribution, the \$k = 1 \ldots n\$th highest out of \$n\$ scores will have a mean equal to

$$\frac{k}{n + 1}$$

i.e. evenly spaced. As \$n\$ increases, the most extreme values \$k = 1, n\$ will get closer and closer to the ends.

For a non-uniform distribution, the distribution of each of these order statistics will be "stretched" according to the inverse CDF. For a "bell-curve" shape like this one, this will tend to stretch out the ends more than the middle. While this stretching doesn't commute with taking the mean, and in RPGs we are typically dealing with discrete rather than continuous distributions, passing the mean through the inverse CDF often gives a decent estimate.

Algorithm notes

While AnyDice is not open-source, the API, performance characteristics, and developer statements imply it is based on enumerating all possible multisets that could come out of a pool. Unfortunately, for large dice, the number of possible multisets can grow quite quickly---in this case, \$\Theta \left( n^{12} \right)\$ where \$n\$ is the number of scores.

Icepool's algorithm is based on a generalization of the technique used in this answer by Ilmari Karonen. It combines the decomposition of a multinomial coefficients into a product of binomials, dynamic programming, and a few other tricks to compute the solution to many types of dice pool problems in polynomial time of reasonable order. In this case, the running time is \$\Theta \left( n^2 \right)\$, times an extra cost of \$\text{O} \left( n^{1.58} \right)\$ if we are seeking an exact fraction.

If you would like to know more, you can read my paper on the subject:

@inproceedings{liu2022icepool,
    title={Icepool: Efficient Computation of Dice Pool Probabilities},
    author={Albert Julius Liu},
    booktitle={Eighteenth AAAI Conference on Artificial Intelligence and Interactive Digital Entertainment},
    volume={18},
    number={1},
    pages={258-265},
    year={2022},
    month={Oct.},
    eventdate={2022-10-24/2022-10-28},
    venue={Pomona, California},
    url={https://ojs.aaai.org/index.php/AIIDE/article/view/21971},
    doi={10.1609/aiide.v18i1.21971}
}

Answer 2

@HighDiceRoller's answer presents what you need (with a nice analysis). dyce¹ is also capable of relatively efficiently² computing the desired results and working around the discussed computation limitations of AnyDice:

from dyce import H, P
d_3_weighted = 3 @ H((2, 3, 3, 4, 5, 6))
p_12_of_3_weighted = 12@P(d_3_weighted)
results_by_selection_index = {
    i: p_12_of_3_weighted.h(i)
    for i in range(0, len(p_12_of_3_weighted))
}

You can play around with a more generalized version in your browser: [source]

While a matter of taste, I find anydyce's³ "burst" graphs often help give a "feel" for distributions like these. (Screenshots below. Note that these appear slightly different from @HighDiceRoller's graphs because anydyce culls data points smaller than Fraction(1, 2**13) for plots. The intention is to increase both readability and performance. The underlying data is the same, however.)

---

¹ dyce is my Python dice probability library.

² In extreme cases that go well beyond what you're trying to do here, like dealing with dozens or hundreds of dice, Icepool is likely to be (often substantially) faster.

³ anydyce is my visualization layer for dyce meant as a rough stand-in for AnyDice.

---

Answer 3

What your AnyDice code calculates is almost certainly not what you're actually rolling at your table.

Basically, your code rolls 9 sets of 3dT (where T is a relabelled d6) and outputs (the distribution of) the highest sum of these 9 sets. Then it rolls another 9 sets of 3dT and outputs (the distribution of) the second highest sum of those, and then it rolls yet another 9 sets of 3dT and outputs (the distribution of) the third highest sum of those, and so on.

I assume that your actual mechanic involves just rolling a single array of 9 (or 12, or however many) sets of 3dT and sorting them.

Now, the AnyDice code you posted does correctly calculate the marginal distributions of the highest / second highest / third highest / etc. stats rolled using this method, but since it rerolls the whole stat array every time, it ends up completely neglecting the dependencies between the stats, such as the fact that the highest stat cannot be lower than the second highest, even though their distributions overlap.

Of course, those dependencies also cannot be represented using a simple set of bar / line charts like AnyDice outputs. Indeed, the full joint distribution that encodes all the dependencies is 9 (or 12, etc.) dimensional and not easily visualizable.

---

Anyway, assuming you want to use AnyDice in "roller mode" to simulate rolling a random array of 3dT stats, all you need is this code:

T: {2,3,3,4,5,6}

output 3dT named "stats"

Then go to the "Roller" tab (the link above will take you there directly), enter 9 (or 12, or any number you want) in the box next to the "Roll" button* and click the button.

The resulting array will not be sorted, so you'll have to sort it yourself. But you were presumably going to manually allocate those values to different stats anyway.

(Using AnyDice for this is honestly kind of overkill, but it does seem to be one of the few free online dice rollers that support custom relabelled dice.)

---

^{*) I don't know how to make an AnyDice link that would prefill the amount box in roller mode with a particular value, and none of the "obvious" URL variations I tried worked. If you know how to do that, please let me know!}

Answer 4

I wrote a little program to solve this using Monte Carlo methods, ie to roll lots of dice and see what the average result is.

Ten million trials seems to be enough to get a results which agree to two decimal places, and a single iteration of that takes about five seconds on my machine. I get this array for twelve stats rolled on 3DT:

7.96, 8.89, 9.58, 10.16, 10.69, 11.19, 11.68, 12.19, 12.75, 13.38, 14.17, 15.36

I know you already have exact answers using various tools, but i thought it would be interesting to see how complicated it is using good old brute force.

Answer 5

I don't think anydice can handle it, but it's not all that hard to calculate instantly in Google Sheets. First I did a quick sheet here to calculate all possibilities of an attribute.

Sum	Count	Prob
6	1	0.46%
7	6	2.78%
8	15	6.94%
9	23	10.65%
10	30	13.89%
11	36	16.67%
12	34	15.74%
13	27	12.50%
14	21	9.72%
15	13	6.02%
16	6	2.78%
17	3	1.39%
18	1	0.46%

From there, we can consider the above table to be one big 216-sided "ABILITY" dice, and then calculate: If we roll 12 ABILITY dice, what are the odds that the Nth largest is exactly M?

The equation for this is the probability that there is 1 die equal to M, anywhere, and 12-N die less than M, and N-1 die greater than M, and the bigger die can be in any of the 11 remaining places. That equation is: P(x=M) * 12 * P(x<M)^(12-N) * P(x>=M)^(N-1) * COMBIN(11, M-1). And I made a Google sheet calculate that:

M\N	1	2	3	4	5	6	7	8	9	10	11	12
6	0	0	0	0	0	0	0	0	0	0	0	0.0556
7	0	0	0	0	0	0	0	0	0	0.0004	0.0162	0.3167
8	0	0	0	0	0	0	0	0.0002	0.0036	0.0358	0.2137	0.58
9	0	0	0	0	0	0.0004	0.0034	0.0214	0.0943	0.2773	0.489	0.392
10	0	0	0	0	0.0037	0.0196	0.0744	0.2019	0.3837	0.486	0.3693	0.1276
11	0	0.0004	0.0034	0.0034	0.0729	0.1919	0.3608	0.4845	0.4555	0.2854	0.1073	0.0183
12	0.0012	0.013	0.0614	0.0614	0.3294	0.4363	0.4127	0.2788	0.1319	0.0416	0.0079	0.0007
13	0.0187	0.1008	0.2468	0.2468	0.355	0.2434	0.1192	0.0417	0.0102	0.0017	0.0002	0
14	0.0952	0.268	0.3428	0.3428	0.1346	0.0482	0.0123	0.0023	0.0003	0	0	0
15	0.2093	0.2744	0.1635	0.1635	0.0139	0.0023	0.0003	0	0	0	0	0
16	0.1979	0.1057	0.0256	0.0256	0.0004	0	0	0	0	0	0	0
17	0.1357	0.0282	0.0027	0.0027	0	0	0	0	0	0	0	0
18	0.0528	0.0027	0.0001	0.0001	0	0	0	0	0	0	0	0
19	0	0	0	0	0	0	0	0	0	0	0	0
20	0	0	0	0	0	0	0	0	0	0	0	0
21	0	0	0	0	0	0	0	0	0	0	0	0

Note: These do not sum to 100%, so I clearly have an error somewhere in my math. You'd think that you could just pick the row of the biggest number in each column as the standard array for that position, but that actually gets subtly wrong answers. But we can use the above table to calculate what we really need: If we roll 12 ABILITY dice, what are the odds that the Nth largest is at least M?

This is a trivial sum, in the same google sheet.

M\N	1	2	3	3	5	6	7	8	9	10	11	12
6	0	0	0	0	0	0	0	0	0	0	0	0.0556
7	0	0	0	0	0	0	0	0	0	0.0004	0.0162	0.3723
8	0	0	0	0	0	0	0	0.0002	0.0036	0.0362	0.2299	0.9523
9	0	0	0	0	0	0.0004	0.0034	0.0216	0.0979	0.3134	0.7189	1.3443
10	0	0	0	0	0.0037	0.02	0.0778	0.2236	0.4816	0.7994	1.0882	1.4719
11	0	0.0004	0.0035	0.0035	0.0766	0.2119	0.4386	0.7081	0.9371	1.0848	1.1955	1.4902
12	0.0013	0.0133	0.0648	0.0648	0.4061	0.6481	0.8513	0.9869	1.0689	1.1264	1.2034	1.4909
13	0.02	0.1141	0.3116	0.3116	0.761	0.8915	0.9705	1.0286	1.0791	1.1281	1.2036	1.4909
14	0.1152	0.3821	0.6544	0.6544	0.8956	0.9397	0.9828	1.0309	1.0794	1.1281	1.2036	1.4909
15	0.3245	0.6565	0.8179	0.8179	0.9096	0.942	0.9831	1.0309	1.0794	1.1281	1.2036	1.4909
16	0.5224	0.7622	0.8435	0.8435	0.9099	0.9421	0.9831	1.0309	1.0794	1.1281	1.2036	1.4909
17	0.6581	0.7903	0.8462	0.8462	0.9099	0.9421	0.9831	1.0309	1.0794	1.1281	1.2036	1.4909
18	0.7109	0.793	0.8462	0.8462	0.9099	0.9421	0.9831	1.0309	1.0794	1.1281	1.2036	1.4909
19	0.7109	0.793	0.8462	0.8462	0.9099	0.9421	0.9831	1.0309	1.0794	1.1281	1.2036	1.4909
20	0.7109	0.793	0.8462	0.8462	0.9099	0.9421	0.9831	1.0309	1.0794	1.1281	1.2036	1.4909
21	0.7109	0.793	0.8462	0.8462	0.9099	0.9421	0.9831	1.0309	1.0794	1.1281	1.2036	1.4909

Then, to get our standard array, we merely pick the row with the value closest to 50% in each column, which I bolded above. Here they are clearly:

16, 15, 14, 13, 12, 12, 11, 11, 10, 9, 9, 7

Again Note: The bottom rows are not 100%, so I clearly have an error somewhere in my math. If anyone mathy can teach me my mistake, that'd be great.

Note#2: My Google Chrome says that sheet is using ~232MB of memory, and calculates ~instantly on my laptop. Far less than you require :P