[OT] Generating distribution of N dice rolls

[Timon Gehr]

On 10.11.22 03:10, H. S. Teoh wrote:
> This is technically OT, but I thought I'd pick the smart brains here for
> my project, which happens to be written in D. ;-)
> 
> Basically, I want to write a function that takes 2 uint arguments k and
> N, and simulates rolling N k-sided dice and counting how many 1's, 2's,
> 3's, ... k's were rolled. Something like this:
> 
> 	uint[k] diceDistrib(uint k)(uint N)
> 		in(k > 0)
> 		in(N > 0)
> 		out(r; r[].sum == N)
> 	{
> 		uint[k] result;
> 		foreach (i; 0 .. N) {
> 			result[uniform(0, k)]++;
> 		}
> 		return result;
> 	}
> 
> The above code works and does what I want, but since N may be large, I'd
> like to refactor the code to loop over k instead of N. I.e., instead of
> actually rolling N dice and tallying the results, the function would
> generate the elements of the output array directly, such that the
> distribution of the array elements follow the same probabilities as the
> above code.
> 
> Note that in all cases, the output array must sum to N; it is not enough
> to merely simulate the roll distribution probabilistically.
> 
> Any ideas?  (Or links if this is a well-studied problem with a known
> solution.)
> ...

The distribution of the output is multinomial(n,k,[1/k,1/k,...,1/k]), 
i.e., n trials, k events, and the events follow an uniform distribution.
https://en.wikipedia.org/wiki/Multinomial_distribution

It is at least as hard to sample as binomial(n,1/k), because that is the 
marginal distribution of each component of the result.

I guess one approach is this (if you can find a way to sample from 
binomial distributions that is good enough for your use case):

auto tot=n;
foreach(i;0..k){
     result[i]=binomial(tot,1.0/(k-i));
     tot-=result[i];
}
return result;

This samples the result one component at a time.

> <ObDReference> I love how D's new contract syntax makes it so conducive
> to expressing programming problem requirements. ;-) </ObDReference>
> 

:)