Setup

import math
from math import pi as pi
from math import cos as cos
from math import sin as sin
import matplotlib
from matplotlib.animation import FuncAnimation as FA
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns

Constants and Helpers

#For animation.
dice_colors = {4: '#FDE725FF',
               6: '#95D840FF',
               8: '#3CBB75FF',
               10: '#33638DFF',
               12: '#453781FF',
               20: '#440154FF'}



#Create dictionary of all dice, name: [numbers]
def getDice() :
    dice = {}
    for f in [4,6,8,10,12,20] :
        dice[f] = [g+1 for g in range(f)]
    return(dice)



dice = getDice()

#Create all possible combinations of four dice.    
def getCombos(dice):
    combos = [[]]
    rolls = 4
    while rolls > 0 :
        new_combos = []
        #Each existing di(c)e combination.
        for c in combos :
            #A new die to that combination.
            for d in dice.keys():
                e = c.copy()
                e.append(d)
                e.sort()
                #Order does not matter.
                if e not in new_combos:
                    new_combos.append(e)
        #Update combos.
        combos = new_combos
        #One less die.
        rolls -= 1
    #Sort least to greatest, by sum, then remove sum.    
    for i in range(126):
        combos[i] = [sum(combos[i])] + combos[i]
    combos.sort()
    for i in range(126):
        combos[i] = combos[i][1:]
    return(combos)



combos = getCombos(dice)

#Frequency distribution from list.
def getFreq(totals) :
    freq = {}
    for t in totals:
        if t in freq.keys() :
            freq[t] += 1
        else:    
            freq[t] = 1
    return(freq)



#Reverse Cumulative frequency distribution from frequency distribution.
def getCumFreq(freq):
    l = len(freq)
    total = sum(freq.values())
    t = sum(freq.values())
    cum_freq = {}
    
    cum_freq[list(freq.keys())[0]] = round(t/total, 7)
    for i in range(1,l) :
        #No longer include previous, but everything else greater.
        t -= freq[list(freq.keys())[i-1]]
        cum_freq[list(freq.keys())[i]] = round(t/total, 7)
    return(cum_freq)



#True dice names.
def diceNames(c) :
    return('d' + str(c[0]) + ', d' + str(c[1]) + ', d' + str(c[2]) + ', and d' + str(c[3]))



#Final output statement.
def printOut(t, d, p) :
    return('Rolling at least a ' + str(t) + ' with ' + diceNames(d) + ' is ' + str(p) + '.')



#Final Outputs.
#If 'fiddler', output the cumulative frequency distribution of 3 d6 and 1 d4.
#If 'ec', output the value and dice combination if the probability is exactly 2/9.
#If 'animation', grab cumulative distribution frequency dictionary for all 126 dice combinations.
def getPossibilities(combos, option) :
    CF = {}
    #x = 0
    for i, c in enumerate(combos) :
        totals = [0]
        #Each dice in combo.
        for d in c :
            new_totals = []
            #Each face on die.
            for f in dice[d] :
                #Each total from previous di(c)e.
                for t in totals:
                    new_totals.append(f+t)
            totals = new_totals
        if option == 'general' :
            pass
            #print(c, len(set(totals)))
        else :
            freq = getFreq(totals)
            cum_freq = getCumFreq(freq)
            CF[i] = cum_freq
            if option == 'fiddler' :
                for cf in cum_freq :
                    print(printOut(cf, combos[0], cum_freq[cf]))
        
            elif option == 'ec' :
                for k in cum_freq.keys() :
                    if cum_freq[k] == 0.2222222:
                        print(printOut(k, c, '2/9'))
    if option == 'animation':
        return(CF)
    

    
#For animation.
#ax - axis.
#v - central vertical position of die.
#d - faces on die.
def getDiceShape(ax, v, d) :
    
    
    #Points along perimeter of die.
    def getPerimeter(o, ps) :
        
        xs = []
        ys = []
        
        #Create perimeter.
        for p in range(ps):
            t = 2*pi*p/ps
            
            #d4 and d6 are larger.
            if d in [4,6] :
                xs.append(5.4*cos(t+o) + 72)
                ys.append((9*1.02/77)*sin(t+o) + v)
                    
            else :
                xs.append(4.8*cos(t+o) + 72)
                #d10 is a squashed hexagon.
                if ((d == 10) and (p not in [0,3])) :
                    ys.append((4*1.02/77)*sin(t+o) + v)
                else :
                    ys.append((8*1.02/77)*sin(t+o) + v)
                
        return(xs+[xs[0]], ys+[ys[0]])

    #Get the perimeter.
    #Also create internal segments.
    #d4 - 3 points.
    if d == 4:
        x, y = getPerimeter(pi/2, 3)

    #d6 - 4 points.
    elif d == 6:
        x, y = getPerimeter(pi/4, 4)

    #d8 - 6 points.   
    elif d == 8: 
        x, y = getPerimeter(pi/2, 6)
        #d8 interior.
        ax.plot([x[0], x[2], x[4], x[0]],
                [y[0], y[2], y[4], y[0]],
                lw = 3,
                c = dice_colors[d])

    #d10 - 6 points (squashed).    
    elif d == 10: 
        x, y = getPerimeter(pi/2, 6)
        #d10 interior.
        ax.plot([x[0], 2/3*x[2]+24, x[2], 2/3*x[2]+24,
                 2/3*x[3]+24, x[3], 2/3*x[3]+24, 
                 2/3*x[4]+24, x[4], 2/3*x[4]+24, x[0]], 
                 [y[0], 2/3*y[2]+v/3, y[2], 2/3*y[2]+v/3,
                 2/3*y[3]+v/3, y[3], 2/3*y[3]+v/3, 
                 2/3*y[4]+v/3, y[4], 2/3*y[4]+v/3, y[0]],
                lw = 3,
                c = dice_colors[d])  
        
    #d12 - 10 points.
    elif d == 12: 
        x, y = getPerimeter(pi/2, 10)
        #d12 interior.
        ax.plot([x[0], 2/3*x[0]+24, 2/3*x[2]+24, x[2],
                 x[2], 2/3*x[2]+24, 2/3*x[4]+24, x[4],
                 x[4], 2/3*x[4]+24, 2/3*x[6]+24, x[6],
                 x[6], 2/3*x[6]+24, 2/3*x[8]+24, x[8],
                 x[8], 2/3*x[8]+24, 2/3*x[0]+24, x[0]],
                [y[0], 2/3*y[0]+v/3, 2/3*y[2]+v/3, y[2],
                 y[2], 2/3*y[2]+v/3, 2/3*y[4]+v/3, y[4],
                 y[4], 2/3*y[4]+v/3, 2/3*y[6]+v/3, y[6],
                 y[6], 2/3*y[6]+v/3, 2/3*y[8]+v/3, y[8],
                 y[8], 2/3*y[8]+v/3, 2/3*y[0]+v/3, y[0]],
                lw = 3,
                c = dice_colors[d])
    
    #d20 - 6 points.
    elif d == 20:
        x, y = getPerimeter(pi/2, 6)
        #d20 interior.
        ax.plot([x[0], 3/4*x[0]+18, x[1], 3/4*x[2]+18,
                 x[2], 3/4*x[2]+18, x[3], 3/4*x[4]+18,
                 x[4], 3/4*x[4]+18, x[5], 3/4*x[0]+18,
                 3/4*x[2]+18, 3/4*x[4]+18, 3/4*x[0]+18],
                [y[0], 3/4*y[0]+v/4, y[1], 3/4*y[2]+v/4,
                 y[2], 3/4*y[2]+v/4, y[3], 3/4*y[4]+v/4,
                 y[4], 3/4*y[4]+v/4, y[5], 3/4*y[0]+v/4,
                 3/4*y[2]+v/4, 3/4*y[4]+v/4, 3/4*y[0]+v/4],
                lw = 3,
                c = dice_colors[d])
        
    #d_ perimeter.
    return(ax.plot(x, y, lw = 4.5, c = dice_colors[d]))



#For animation.
#ax - axis.
#v - vertical position.
#text - text.
#s - size.
#w - weight.
def annotationText(ax, v, text, s, w) :
    return(ax.text(54, v,
                   text,
                   ha = 'center',
                   va = 'top',
                   size = s,
                   weight = w,
                   c = 'k'))



#For animation.
#ax - axis.
#v - central vertical position of die.
#d - faces on die.
def getDiceText(ax, v, d):
    return(ax.text(72, v, d,
                   ha = 'center',
                   va = 'center',
                   size = 25,
                   weight = 'bold',
                   c = dice_colors[d]))

Answers

getPossibilities([combos[6]], 'fiddler')

Rolling at least a 4 with d4, d6, d6, and d6 is 1.0.
Rolling at least a 5 with d4, d6, d6, and d6 is 0.9988426.
Rolling at least a 6 with d4, d6, d6, and d6 is 0.994213.
Rolling at least a 7 with d4, d6, d6, and d6 is 0.9826389.
Rolling at least a 8 with d4, d6, d6, and d6 is 0.9594907.
Rolling at least a 9 with d4, d6, d6, and d6 is 0.9201389.
Rolling at least a 10 with d4, d6, d6, and d6 is 0.8599537.
Rolling at least a 11 with d4, d6, d6, and d6 is 0.7777778.
Rolling at least a 12 with d4, d6, d6, and d6 is 0.6759259.
Rolling at least a 13 with d4, d6, d6, and d6 is 0.5601852.
Rolling at least a 14 with d4, d6, d6, and d6 is 0.4398148.
Rolling at least a 15 with d4, d6, d6, and d6 is 0.3240741.
Rolling at least a 16 with d4, d6, d6, and d6 is 0.2222222.
Rolling at least a 17 with d4, d6, d6, and d6 is 0.1400463.
Rolling at least a 18 with d4, d6, d6, and d6 is 0.0798611.
Rolling at least a 19 with d4, d6, d6, and d6 is 0.0405093.
Rolling at least a 20 with d4, d6, d6, and d6 is 0.0173611.
Rolling at least a 21 with d4, d6, d6, and d6 is 0.005787.
Rolling at least a 22 with d4, d6, d6, and d6 is 0.0011574.

getPossibilities(combos, 'ec')

Rolling at least a 16 with d4, d6, d6, and d6 is 2/9.
Rolling at least a 21 with d6, d6, d8, and d10 is 2/9.
Rolling at least a 36 with d10, d12, d12, and d20 is 2/9.

Animation

CF = getPossibilities(combos, 'animation')

sns.set()
fig = plt.figure(figsize = (15, 9))

#Animate.
def animate(i):
    
    ax = fig.add_subplot(xlim = (3.5, 80.5),
                         ylim = (0, 1.02))

    #Title setup.
    ax.set_title('Cumulative Probability of Greater Than or Equal to Four Dice Sum', fontsize = 24)

    #X-axis setup.
    ax.set_xlabel('Sum', fontsize = 20)
    ax.set_xticks([4*x for x in range(1, 21)])

    #Y-axis setup.
    ax.set_ylabel('Cumulative Probability', fontsize = 20)
    ax.tick_params(axis = 'both', which = 'major', labelsize = 16)
    
    #Combination.
    annotationText(ax, 0.98, "Combination :", 24, 'normal')
    
    #Dice Sum.
    annotationText(ax, 0.78, "Max Dice Sum :", 24, 'normal')
    
    #Data from specific dice combination.
    data = CF[i]
    k = list(data.keys())
    v = list(data.values())
    
    #All bars.
    bars0 = plt.bar(x = k,
                    height = v,
                    color = '#21908CFF',
                    edgecolor = 'k')
    
    #2/9.
    if 0.2222222 in v:
        j = v.index(0.2222222)
        barsA = plt.bar(x = k[j],
                        height = v[j],
                        color = 'k',
                        edgecolor = 'k')
        textA = ax.text(k[j]+2, v[j],
                        '2/9',
                        ha = 'left',
                        va = 'bottom',
                        size = 30,
                        weight = 'bold',
                        c = 'k')
        
    
    #Die 0.
    shape0 = getDiceShape(ax, 0.86, combos[i][0])
    text0 = getDiceText(ax, 0.855, combos[i][0])

    #Die 1.
    shape1 = getDiceShape(ax, 0.63, combos[i][1])
    text1 = getDiceText(ax, 0.625, combos[i][1])

    #Die 2.
    shape2 = getDiceShape(ax, 0.4, combos[i][2])
    text2 = getDiceText(ax, 0.395, combos[i][2])

    #Die 3.
    shape3 = getDiceShape(ax, 0.17, combos[i][3])
    text3 = getDiceText(ax, 0.165, combos[i][3])
    
    #Combo.
    text4 = annotationText(ax, 0.92, i+1, 30, 'bold')
    
    #Sum.
    text5 = annotationText(ax, 0.72, sum(combos[i]), 30, 'bold')
    
    return bars0, shape0, text0, shape1, text1, shape2, text2, shape3, text3, text4, text5

#Run animation.
anim = FA(fig, animate, frames = 126, interval = 400)   

#Save animation.
anim.save('2024.11.29.mp4');

Appendix

dice

{4: [1, 2, 3, 4],
 6: [1, 2, 3, 4, 5, 6],
 8: [1, 2, 3, 4, 5, 6, 7, 8],
 10: [1, 2, 3, 4, 5, 6, 7, 8, 9, 10],
 12: [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12],
 20: [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20]}

for i, c in enumerate(combos):
    print(str(i) + ' : ' + str(c))

0 : [4, 4, 4, 4]
1 : [4, 4, 4, 6]
2 : [4, 4, 4, 8]
3 : [4, 4, 6, 6]
4 : [4, 4, 4, 10]
5 : [4, 4, 6, 8]
6 : [4, 6, 6, 6]
7 : [4, 4, 4, 12]
8 : [4, 4, 6, 10]
9 : [4, 4, 8, 8]
10 : [4, 6, 6, 8]
11 : [6, 6, 6, 6]
12 : [4, 4, 6, 12]
13 : [4, 4, 8, 10]
14 : [4, 6, 6, 10]
15 : [4, 6, 8, 8]
16 : [6, 6, 6, 8]
17 : [4, 4, 8, 12]
18 : [4, 4, 10, 10]
19 : [4, 6, 6, 12]
20 : [4, 6, 8, 10]
21 : [4, 8, 8, 8]
22 : [6, 6, 6, 10]
23 : [6, 6, 8, 8]
24 : [4, 4, 10, 12]
25 : [4, 6, 8, 12]
26 : [4, 6, 10, 10]
27 : [4, 8, 8, 10]
28 : [6, 6, 6, 12]
29 : [6, 6, 8, 10]
30 : [6, 8, 8, 8]
31 : [4, 4, 4, 20]
32 : [4, 4, 12, 12]
33 : [4, 6, 10, 12]
34 : [4, 8, 8, 12]
35 : [4, 8, 10, 10]
36 : [6, 6, 8, 12]
37 : [6, 6, 10, 10]
38 : [6, 8, 8, 10]
39 : [8, 8, 8, 8]
40 : [4, 4, 6, 20]
41 : [4, 6, 12, 12]
42 : [4, 8, 10, 12]
43 : [4, 10, 10, 10]
44 : [6, 6, 10, 12]
45 : [6, 8, 8, 12]
46 : [6, 8, 10, 10]
47 : [8, 8, 8, 10]
48 : [4, 4, 8, 20]
49 : [4, 6, 6, 20]
50 : [4, 8, 12, 12]
51 : [4, 10, 10, 12]
52 : [6, 6, 12, 12]
53 : [6, 8, 10, 12]
54 : [6, 10, 10, 10]
55 : [8, 8, 8, 12]
56 : [8, 8, 10, 10]
57 : [4, 4, 10, 20]
58 : [4, 6, 8, 20]
59 : [4, 10, 12, 12]
60 : [6, 6, 6, 20]
61 : [6, 8, 12, 12]
62 : [6, 10, 10, 12]
63 : [8, 8, 10, 12]
64 : [8, 10, 10, 10]
65 : [4, 4, 12, 20]
66 : [4, 6, 10, 20]
67 : [4, 8, 8, 20]
68 : [4, 12, 12, 12]
69 : [6, 6, 8, 20]
70 : [6, 10, 12, 12]
71 : [8, 8, 12, 12]
72 : [8, 10, 10, 12]
73 : [10, 10, 10, 10]
74 : [4, 6, 12, 20]
75 : [4, 8, 10, 20]
76 : [6, 6, 10, 20]
77 : [6, 8, 8, 20]
78 : [6, 12, 12, 12]
79 : [8, 10, 12, 12]
80 : [10, 10, 10, 12]
81 : [4, 8, 12, 20]
82 : [4, 10, 10, 20]
83 : [6, 6, 12, 20]
84 : [6, 8, 10, 20]
85 : [8, 8, 8, 20]
86 : [8, 12, 12, 12]
87 : [10, 10, 12, 12]
88 : [4, 10, 12, 20]
89 : [6, 8, 12, 20]
90 : [6, 10, 10, 20]
91 : [8, 8, 10, 20]
92 : [10, 12, 12, 12]
93 : [4, 4, 20, 20]
94 : [4, 12, 12, 20]
95 : [6, 10, 12, 20]
96 : [8, 8, 12, 20]
97 : [8, 10, 10, 20]
98 : [12, 12, 12, 12]
99 : [4, 6, 20, 20]
100 : [6, 12, 12, 20]
101 : [8, 10, 12, 20]
102 : [10, 10, 10, 20]
103 : [4, 8, 20, 20]
104 : [6, 6, 20, 20]
105 : [8, 12, 12, 20]
106 : [10, 10, 12, 20]
107 : [4, 10, 20, 20]
108 : [6, 8, 20, 20]
109 : [10, 12, 12, 20]
110 : [4, 12, 20, 20]
111 : [6, 10, 20, 20]
112 : [8, 8, 20, 20]
113 : [12, 12, 12, 20]
114 : [6, 12, 20, 20]
115 : [8, 10, 20, 20]
116 : [8, 12, 20, 20]
117 : [10, 10, 20, 20]
118 : [10, 12, 20, 20]
119 : [4, 20, 20, 20]
120 : [12, 12, 20, 20]
121 : [6, 20, 20, 20]
122 : [8, 20, 20, 20]
123 : [10, 20, 20, 20]
124 : [12, 20, 20, 20]
125 : [20, 20, 20, 20]

Rohan Lewis

2024.12.02