Setup

import math
import matplotlib.pyplot as plt
from math import factorial as MF
import numpy as np
import pandas as pd
import seaborn as sns

prob_dist = [1/8, 3/8, 3/8, 1/8]

abcds = [(8, 0, 0, 0),
         (7, 1, 0, 0),
         (6, 2, 0, 0),
         (6, 1, 1, 0),
         (5, 3, 0, 0),
         (5, 2, 1, 0),
         (5, 1, 1, 1),
         (4, 4, 0, 0),
         (4, 3, 1, 0),
         (4, 2, 2, 0),
         (4, 2, 1, 1),
         (3, 3, 2, 0),
         (3, 3, 1, 1),
         (3, 2, 2, 1),
         (2, 2, 2, 2)]

def overcounting(x) :
    freq = {}
    for i in x :
        if i not in freq.keys():
            freq[i] = 1
        else :
            freq[i] += 1
    
    oc = 1
    
    for f in freq.values():
        oc *= MF(f)
    
    return(oc)

probs = []
count = 0
for abcd in abcds :
    
    a = abcd[0] 
    b = abcd[1] 
    c = abcd[2] 
    d = abcd[3]
    prob = 0
    
    for w in prob_dist :
        pd1 = prob_dist.copy()
        pd1.remove(w) 
        for x in pd1 :
            pd2 = pd1.copy()
            pd2.remove(x)
            for y in pd2 :
                pd3 = pd2.copy()
                pd3.remove(y)
                z = pd3[0]
                
                prob += (w**a) * (x**b) * (y**c) * (z**d) * MF(8) / (MF(a) * MF(b) * MF(c) * MF(d))
                
    probs.append(round(prob / overcounting(abcd), 5))
    print(abcd, overcounting(abcd))
    
    
df = pd.DataFrame({"X" : range(1,16),
                   "Quadruple": abcds,
                   "Probability": probs})

df

(8, 0, 0, 0) 6
(7, 1, 0, 0) 2
(6, 2, 0, 0) 2
(6, 1, 1, 0) 2
(5, 3, 0, 0) 2
(5, 2, 1, 0) 1
(5, 1, 1, 1) 6
(4, 4, 0, 0) 4
(4, 3, 1, 0) 1
(4, 2, 2, 0) 2
(4, 2, 1, 1) 2
(3, 3, 2, 0) 2
(3, 3, 1, 1) 4
(3, 2, 2, 1) 2
(2, 2, 2, 2) 24

	X	Quadruple	Probability
0	1	(8, 0, 0, 0)	0.00078
1	2	(7, 1, 0, 0)	0.01044
2	3	(6, 2, 0, 0)	0.02683
3	4	(6, 1, 1, 0)	0.03417
4	5	(5, 3, 0, 0)	0.04741
5	6	(5, 2, 1, 0)	0.12809
6	7	(5, 1, 1, 1)	0.02956
7	8	(4, 4, 0, 0)	0.02873
8	9	(4, 3, 1, 0)	0.17504
9	10	(4, 2, 2, 0)	0.08201
10	11	(4, 2, 1, 1)	0.12797
11	12	(3, 3, 2, 0)	0.08532
12	13	(3, 3, 1, 1)	0.07090
13	14	(3, 2, 2, 1)	0.14059
14	15	(2, 2, 2, 2)	0.01217

df["Probability"].sum()

1.0000099999999998

sns.set()
fig = plt.figure(figsize = (16, 10))
ax = fig.add_subplot(xlim = (0.5, 15.5),
                     ylim = (0, 0.2))


#Title setup.
ax.set_title('Probabilities of Distributions of 3 Coin Flips 8 Times', fontsize = 24)

#X-axis setup.
ax.set_xlabel('(a, b, c, d)', fontsize = 22)
ax.set_xticks(range(1, 16))
ax.set_xticklabels(abcds, rotation = 65, ha = "right")
#Y-axis setup.
ax.set_ylabel('Probability', fontsize = 22)
ax.tick_params(axis = 'both', which = 'major', labelsize = 18)


bars0 = plt.bar(x = df['X'],
                height = df['Probability'],
                color = '#21908CFF',
                width = 0.7,
                edgecolor = 'k')

#Labels.
for i in range(15):
    plt.text(x = bars0[i].get_xy()[0] + 0.33,
             y = bars0[i].get_height() + .002,
             s = bars0[i].get_height(),
             ha = 'center',
             size = 'large',
             weight = 'bold',
             c = '#481567FF')
    
    
#Save.
fig.savefig("2023.11.17 Extra Credit.png",
            bbox_inches = 'tight')

Rohan Lewis

2023.11.17