XKCD 3015

xkcd.py

from sympy import *
import functools

a, b, z = var('a b z')

# discrete uniform pgf
G_du = (z ** a - z ** (b + 1)) / ((b - a + 1) * (1 - z))

# pgf for 3d6 + 1d4
G = G_du.subs(a, 1).subs(b, 6) ** 3 * G_du.subs(a, 1).subs(b, 4)

@functools.cache
def Gk(k):
    """k-th derivative of G"""
    if k == 0:
        return simplify(G)
    elif k > 0:
        return simplify(diff(Gk(k - 1), z))

def p(k):
    """pmf of G"""
    return Gk(k).subs(z, 0) / factorial(k)

for k in range(22 + 1):
    cmf = sum([p(l) for l in range(0, k + 1)])
    print(k, p(k), cmf)
    # note that k=16 has cumulative mass 7/9