Background
I found this interesting question Formula for dropping dice (non-brute force) and excellent answer https://stats.stackexchange.com/a/242857/221422, but couldn't figure out how to generalize a generating function for when more than one die is dropped. Similarly, I'm having difficulty figuring out a related mechanic for when the highest roll is dropped.
Description of the Problem
Suppose you have $N$ fair dice each with $S$ sides. Roll all the dice and then remove the lowest [or highest, alternatively] $M$ (where $M > 0$ and $M < N$) dice and then sum the remainder. What is the probability distribution of the sum? Specifically, how does one go about finding the generating polynomial?
Implementation of whuber's answer
I found whuber's answer to be incredibly thorough. I thought it might be nice to see how to actually implement it in code, so I've pasted that below.
from functools import reduce
from numpy.polynomial import polynomial as p
def generating_function(k, d, n):
return p.polypow(
[0] * k + [1] * (d - k + 1),
n
)
def drop_one_die(n, d):
tmp = [
generating_function(k, d, n) for k in range(1, d + 2)
]
differences = (
(tmp[i] - tmp[i + 1])[i + 1:] for i in range(d)
)
return reduce(p.polyadd, differences)
print(
drop_one_die(4, 6)
)
Other considerations / Multinomial distribution
To generalize even further, instead of a fair die where each outcome is equally likely, what if you start with a general multinomial distribution?
So instead of
$$(1/6)x + (1/6)x^2 + (1/6)x^3 + (1/6)x^4 + (1/6)x^5 + (1/6)x^6$$
you start with
$$p_0 + {p_1}{x} + {p_2}{x^2} + ... + {p_n}{x^n}$$
Thanks!