# 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?

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)
)

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?

$$(1/6)x + (1/6)x^2 + (1/6)x^3 + (1/6)x^4 + (1/6)x^5 + (1/6)x^6$$
$$p_0 + {p_1}{x} + {p_2}{x^2} + ... + {p_n}{x^n}$$